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I.  INTRODUCTION 

One  of  the  most  critical  functions  of  a  submersible  is  accurate  depth  keeping  at  the 
commanded  depth  Such  a  function  can  be  carried  out  either  manually  or  automatically, 
especially  in  cases  where  human  intervention  is  not  possible  or  desirable  Due  to  the 
technological  significance  and  numerous  scientific  applications  of  submersible  vehicle 
systems,  design  issues  of  appropriate  depth  keeping  control  laws  have  received  wide 
attention  in  the  past  Such  control  system  designs  include  linear  and  nonlinear  controllers 
[Refs  1&2],  model  based  compensators  [Ref  3],  adaptive  control  [Ref  4],  and  sliding 
mode  control  laws  [Refs.  5&6] 

Response  accuracy  and  stability  are  primary  considerations  in  designing  depth 
keeping  control  law.  Of  paramount  importance,  in  this  area,  are  the  robustness  properties 
of  the  particular  design;  i.e.,  its  ability  to  maintain  accuracy  and  stability  in  the  presence  of 
incomplete  sensor  and  environmental  information,  as  well  as  actual/mathematical  model 
mismatch  The  scope  of  the  work  in  this  thesis  is  to  demonstrate  a  potential  loss  of 
stability  that  may  occur  when  a  submarine  is  operating  at  low  speeds  This  loss  of  stability 
can  occur  regardless  of  the  particular  means  used  for  depth  control.  The  study  is 
accomplished  through  the  use  of  an  eigenvalue  analysis,  a  steady  state  analysis  and  a 
controllability  analysis  [Ref.  7].  It  is  shown  that  such  a  loss  of  stability  is  accompanied  by  a 
slow  divergence  of  trajectories  away  from  the  commanded  path  Solution  branching 
occurs  in  the  form  of  generic  pitchfork  bifurcations  [Refs  8,  9].  A  complete 
characterization  of  the  problem  is  given  utilizing  singularity  techniques,  which  have  been 
proven  to  be  very  useful  in  the  analysis  of  similar  problems  [Refs.  10,  11,  12].  The  use  of 


bifurcation  theory  allows  the  crucial  vehicle  parameters  that  govern  the  problem  of 
solution  branching  to  be  determined,  and  the  develop  guidelines  to  prevent  its  occurrence 

Finally,  a  new  look  at  the  problem  of  dive  plane  reversal  [Ref  13],  based  on  solution 
branching  results  is  presented  The  term  dive  plane  reversal  refers  to  a  well-known 
phenomenon  in  submarine  operations  where,  during  low  speed  depth  keeping,  there  is  a 
need  to  reverse  the  direction  of  dive  plane  deflection  in  order  to  execute  a  given  change  in 
depth  Physically,  this  can  be  explained  by  considering  the  relative  magnitude  of  the 
hydrodynamic  forces  At  moderate  and  high  speeds,  the  normal  force  on  the  submarine's 
hull  due  to  the  angle  of  attack  exceeds  the  normal  dive  plane  force  and  the  boat  responds 
to  ordered  dive  plane  angles  as  expected.  The  phenomenon  of  dive  plane  reversal  occurs 
at  speeds  below  a  certain  critical  speed  in  which  the  normal  hull  force  is  less  than  the 
normal  dive  plane  force  and  the  response  of  the  boat  is  reversed.  [Ref.  14]  Vehicle 
modeling  in  this  work  follows  standard  notation  [Ref.  15]  and  numerical  results  are 
presented  for  the  DARPA  SUBOFF  model  [Ref  16]  for  which  a  set  of  hydrodynamic 
coefficients  and  geometric  properties  is  available.  Special  emphasis  is  given  to  identifying 
the  proper  non-dimensional  parameters  in  the  problem,  so  that  extension  of  these  results 
to  full  scale  models  and  other  designs  is  possible  using  minimal  experimental  and/or 
analytical  results 


D.  VEHICLE  MODELING 


A.     EQUATIONS  OF  MOTION 

1.  Introduction 

For  the  purpose  of  this  problem,  and  the  subject  of  the  maneuvering  and  motion 
control  of  the  vehicle,  the  following  assumptions  are  made: 

1 )  The  vehicle  behaves  as  a  rigid  body, 

2)  The  earth's  rotation  is  negligible  as  far  as  acceleration  components  of  the 
center  of  mass  are  concerned, 

3)  The  primary  forces  that  act  on  the  vehicle  have  inertial,   gravitational, 
hydrostatic  and  hydrodynamic  origins. 

2.  Coordinate  Systems  And  Positional  Definitions 

A  global  navigation  frame,  OXYZ,  as  shown  in  Figure  2.1,  is  defined  with  origin, 
O,  and  a  set  of  axes  aligned  with  directions  North,  East  and  Down  This  produces  a  right- 
hand  reference  frame  with  unit  vectors  7,  J,  and  K  Ignoring  the  earth's  rotation  rate  in 
comparison  to  the  angular  rates  produced  by  the  vessel's  motion,  it  can  be  said  that  the  I , 
J,  and  K  coordinate  frame  is  an  inertial  reference  frame  in  which  Newton's  Laws  of 
Motion  will  be  valid  As  seen  in  Figure  2  1,  a  vehicle's  position  Ro  ,  in  this  frame  will  have 
the  vector  components,        R    -  [ XJ  +Y0J  +  Z0KJ .  A  standard  convention  that  is  used 

in  both  aircraft  and  marine  vehicle  dynamics  places  the  Y  axis  to  the  right  while  looking 
along  the  X  axis,  and  the  Z  axis  is  positive  downwards  Figure  2.1  also  shows  a  vehicle 
with  some  general  attitude  in  the  original  coordinate  system. 


Figure  2-1.  Coordinate  Axes  Convention. 
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Now  define  a  body  fixed  frame  of  reference  Oxyz ,  with  the  origin  O  ,  located  on 
the  vehicle  centerline,  moving  and  rotating  with  the  vehicle,  in  which  the  vehicle's  center 
of  mass,  G,  has  some  position  other  than  the  origin  of  the  vehicle  fixed  frame  (refer  to 
Figure  2-1)  This  origin  will  be  the  point  about  which  all  vehicle  body  force  will  be 
computed  in  later  sections  of  this  chapter  The  convention  used  for  the  vehicle  fixed 
frame,  with  unit  vectors  /  ,  j  ,  k ,  is  that  its  origin  lies  at  the  ship's  center  (origin  at  the  half 
length),  the  x  axis  at  main  deck  level,  parallel  to  the  longitudinal  centerline,  with  the  z  axis 
vertically  down  The  vessel's  center  of  gravity  and  center  of  buoyancy  do  not  generally  lie 
at  the  origin  of  the  body  fixed  frame,  nor  are  they  collocated.  These  points  are  denoted  by 
G  and  B  respectively.  The  vector  components  of  pG  and  pg  are  thus  [xj  +  yGj  +  zGk  /, 
and  [xgT  +  yBj  +  zBk  ]  .  The  location  of  the  center  of  mass  is  important  because 
Newton's  Laws  of  Motion  equate  forces  on  a  body  to  the  rate  of  change  of  linear 
momentum  of  the  center  of  mass  and  moments  about  the  body  center  of  mass  to  the  rate 
of  change  of  angular  momentum.  This  is  a  particularly  important  point  to  bear  in  mind  for 
ship  and  submarine  motion  dynamics  as  the  center  of  buoyancy  is  determined  by  the  shape 
of  the  submerged  portion  of  the  ship  body  while  the  center  of  gravity  is  determined  by  the 
distribution  of  the  weight  over  the  entire  ship  body.  Having  defined  a  coordinate  system 
that  will  be  used  to  describe  a  ship's  position,  there  is  a  need  to  define  angular  orientation 
of  the  ship's  body,  leading  to  the  definition  of  translational  and  rotational  velocities  and 
accelerations 

3.      Angular  Position  In  The  Global  Reference  Frame 

There  are  several  different  ways  that  the  attitude  of  a  vehicle  can  be  described  in 
reference  to  the  global  frame.  The  most  usual  and  common  method  is  to  define  three 
angles,  called  Euler  angles  that  uniquely  define  the  angular  orientation  of  the  vehicle 
reference  frame,  relative  to  the  global  reference  frame  One  problem  with  the  Euler  angle 


approach  is  that  a  singularity  exists  when  one  of  the  angles  reaches  90  degrees  (an  aircraft 
in  pure  vertical  flight  cannot  distinguish  its  azimuth  angle  from  its  roll  angle).  This 
limitation  which  can  sometimes,  although  rarely,  cause  trouble  in  flight  simulations  and 
control  computations  can  be  overcome  by  the  use  of  quaternions  which  introduce  four 
rather  than  three  variables  to  describe  angular  position.  In  this  presentation,  however,  we 
will  useEuler  angles  as  it  is  the  most  widely  used  method,  although  the  use  of  quaternions 
has  found  favor  in  robotics  applications  and  computer  graphics. 

While  any  consistent  definition  of  three  base  angles  would  be  sufficient,  the  most 
convenient  formulation  of  these  angles  is  in  the  widely  used  military  terms  of  an  azimuth, 
elevation,  and  spin  notation  The  rates  of  change  of  these  angles  do  not  generally 
correspond  to  the  other  commonly  used  angular  rates  describing  angular  velocity  of  a 
vehicle,  that  is  yaw  rate,  pitch  rate  and  roll  rate  except,  as  will  be  seen,  where  motion  is 
limited  to  small  angle  rotations 

4.      Rotational  Transformations 

Vehicle  attitude  is  important  in  defining  position.  Under  dynamic  conditions,  the 
pitch  or  roll  can  cause  problems  for  manned  vessels  and  the  heading  is  critical  in 
navigation.  For  the  purpose  of  considering  angular  rotations,  consider  a  forward 
transformation  from  a  coordinate  triad  aligned  with  the  global  reference  frame  and  perform 
a  sequence  of  three  rotations  to  finally  align  the  result  with  a  frame  that  is  assumed  to  be 
parallel  to  the  vehicle  body  coordinate  axes,  and  moving  with  the  vehicle  at  all  times 
Begin  by  defining  an  azimuth  rotation,  \\i ,  as  a  positive  rotation  about  the  global  Z  axis. 
Next  define  a  subsequent  rotation  6 ,  (positive  up)  about  the  new  Y  axis,  followed  by  a 
positive  rotation  <J> ,  about  the  new  X  axis  The  triple  rotational  transformation  in  terms  of 

these  three  angles,  is  then  sufficient  to  describe  the  angular  orientation  of  the  vehicle  at  any 
time    It  follows  that  any  position  vector,  R0,  in  an  original  reference  frame  given  by 


Ro  =  [x0,y0,z0] ,  will  have  different  coordinates  in  a  rotated  frame  when  an  azimuth 
rotation  by  angle  vj/,  is  made  about  the  global  Z  axis. 

If  the  new  position  is  defined  by,  .ft,  =  [x],y],z]  /,  it  can  be  seen  that  there  is  a 
relation  between  the  vector's  coordinates  in  the  new  reference  frame  with  those  that  it  had 
in  the  old  reference  frame  Using  trigonometrical  relationships  and  Figure  2.2  as  a  guide,  it 
follows  that, 

X]  =  X0  cosy  +  Y0  sin  y  (2.1) 

Yx  =  - X0  sin  y  +  Y0  cosy  (2.2) 

This  relation  can  be  expressed  in  matrix  form  by  the  rotation  matrix  operation, 

%=[T,.Z]R0,  (2-3) 

where  the  rotation  matrix  \TViZ\  represents  an  orthogonal  transformation. 
Premultiplication  of  this  rotation  matrix  with  any  vector,  R0,  will  produce  the  components 
of  the  same  vector  in  the  rotated  coordinate  frame.  Continuing  with  the  series  of  rotations 
results  in  a  combined  total  rotation  transformation, 

T(W,y)  =  m)T(Q)T(y).  (2.4) 

In  expanded  notation  equation  2.4  takes  the  form; 

cos  y  cos  6  sin  y  cos  6  -  sin  6 

cos  \\f  sin  6  sin  <f)  -  sin  \\f  cos  §    sin  \j/  sin  9  sin  §  +  cos  \\f  cos  §    cos  6  sin  (J) 
cosy  sin  6  cos  (J)+ sin  \\f  sin  <f>    siny  sin  Q  cos  §- cosy  sin  §    cos  Q  cos  § 

and  it  can  be  said  that  any  position  vector  in  an  original  reference  frame  may  be 
expressed  in  a  rotated  frame  with  coordinates  given  by  the  operation, 

/L=[W,e,y;Rw.  (2.5) 


xAj  1 


Figure  2-2.  Azimuth  Rotation. 


5.      Kinematics 

Kinematics  deals  with  the  relationships  of  motion  quantities  regardless  of  the 
forces  induced  by  their  prescribed  motions  Descriptions  of  a  ship's  position  both 
translational  and  rotational,  will  need  to  be  related  to  velocities,  both  translational  and 
rotational,  prior  to  extending  the  situation  to  accelerations.  The  connection  between 
translational  velocity  and  the  rate  of  change  of  translational  position  is  straight  forward 
Define  the  global  velocity  as, 


R  = 


X 
Y 
Z 


(2.6) 


This  vector  will  have  components  that  are  different  when  seen  in  a  body  fixed  frame 
Define  the  body  fixed  components  of  the  global  velocity  vector  as  [u,  v,  w]\  These 
components,    in    terms    of  the    above    global    quantities    are    given    by    the    forward 
transformation  defined  earlier  to  be, 


(2.7) 


It  is  a  simple  reverse  coordinate  transformation  from  body  fixed  to  global  coordinates  to 
see  that. 


// 

X 

V 

=  T(<\>,Q,y) 

Y 

w 

Z 

X 
Y 
Z 


=  rY<i>,e,vj/; 


(2.8) 


This  shows  that  the  progression  of  a  vehicle  in  a  global  frame  clearly  depends  on  its  local 
velocity  components  and  its  attitude  Put  simply,  u  corresponding  to  a  vehicle's  forward 


v 


speed  (surge),  v  corresponding  to  a  side  slip  velocity  (sway),  and  w  corresponding  to  any 
velocity  component  in  the  local  Z  direction  (heave)  and  the  vehicle's  global  velocity 
components  depend  on  heading,  pitch,  and  roll  attitude. 

The  connection  between  angular  attitude  and  angular  velocity  is  not  quite  so 
apparent  At  first  sight,  it  is  tempting  to  define  the  instantaneous  angular  velocity  of  the 
vehicle  simply  as  the  rate  of  change  of  its  angular  position  defined  by  the  Euler  angles. 
This  is  erroneous  however,  because  the  rotation  0 ,  was  defined  as  a  rotation  about  the 
intermediate  frame  after  a  rotation  y  had  been  made  Vehicle  inertial  angular  rates  are 
defined  in  terms  of  components  that  have  angular  velocities  about  the  global  axes.  It  is 
necessary  to  relate  both  Euler  angles  and  their  rates  of  change  to  angular  velocity 
components  about  the  global  axes  to  their  components  lying  along  the  body  fixed  axes  in 
any  attitude  The  prime  reason  for  this  is  that  it  is  difficult  to  construct  physical  sensors  to 
measure  rates  of  change  of  Euler  angles  However,  rate  gyros  in  common  use  today  are 
easily  constructed  to  measure  the  components  of  the  inertial  angular  velocity  of  a  vehicle 
that  lie  along  the  vehicle's  body  axes  It  follows  that  the  instantaneous  angular  velocity  of 
the  vehicle  can  be  related  to  the  instantaneous  rate  of  change  of  angular  orientation  only 
after  considerations  of  the  intermediate  transformations  used.  In  other  words,  if  a  is 
defined  as  the  angular  attitude  vector,  a  =  [\|/,6,<J>] ;  and  the  inertial  components  of  the 
vehicle  angular  rate  lying  along  the  body  axes  as  co  =  [p,q,r] ;  then,  d  =  /(($)  The  details 

of  the  nonlinear  functional  relations  involved  are  provided  by  viewing  the  rate  of  change  of 
the  rotation  y  as  a  vector  quantity  lying  along  the  original  Z  axis.  The  rate  of  change  of 
the  angle  6  is  viewed  as  a  vector  quantity  lying  along  the  Y  axis  of  the  first  intermediate 
frame,  and  the  rate  of  change  of  the  angle  <})  is  viewed  as  a  vector  lying  along  the  X  axis  of 
the  final  (body  fixed)  frame.  Each  of  these  component  rates  of  change  of  angular  position 
has  component  parts  that  project  onto  the  final  frame  and  it  is  the  sum  total  of  all  the 
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components  that  give  the  total  angular  velocity  as  seen  in  the  final  frame  of  reference 
Using  the  required  transformations  for  the  rate  components  from  each  Euler  angle,  we  get, 


=  2,T(VT(B)T(y) 


"o" 

V 

+  m)T(Q) 

e 

+  m) 

0 

_o_ 

0 

(2.9) 


with  the  result; 


-vj/s//70  +  <j> 

VJ/  57/70 +  0C0S<t> 

\y  cos  0  cos  <{>  -  0  sin  (j) 


(2.10) 


Inverting  equation  2.10  yields  a  solution  for  the  rates  of  change  of  the  Euler  angles  in 
terms  of  the  body  fixed  components  of  the  angular  velocity  vector, 


p  +  qsin§tanQ  +  r  costy  tanQ 

qcos§-r  sin§ 

(q  sin  <f)  +  r  cos§)  /  cosQ 


(2.11) 


Notice  that  for  small  angular  rotations,  as  expected, 

e  =  q; 
\\f  =  r. 

At  this  point  the  kinematic  relationships  between  velocities,  as  seen  in  the  body  fixed 

frame,  and  the  rates  of  change  of  global  positions  and  Euler  angles  have  been  defined  The 

resulting  set  of  differential  equations  forms  a  consistent  set  in  that  given  a  set  of  vehicle 

velocity  data  versus  time,  its  position  and  attitude  may  be  computed 
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6.  Translational  Equations  Of  Motion 

The  global  acceleration  of  the  center  of  mass  is  derived  by  differentiating  the 
velocity  vector,  R^,  taking  into  account  that  the  center  of  mass  lies  in  a  rotating  reference 
frame  Considering  the  total  differentiation,  the  global  acceleration  of  the  center  of  mass 
becomes, 

^,  =  v  +  ci)xpG+a)XG)xpG  +coxv,  (2. 12) 

where  v  =  Ro  The  translational  equation  of  motion  is  found  by  equating  this  acceleration 
to  the  net  sum  of  all  forces  acting  on  the  vehicle  in  three  degrees  of  freedom  (X,Y,Z).  One 
important  factor  to  recognize  is  that  the  equation  of  motion  derived  in  this  manner  is  a 
vector  equation  with  the  components  expressed  in  the  body  fixed  frame  and  unit  vectors 
I,  j  and  k  This  has  been  deliberately  done  because  the  dominant  forces  acting  on  a 
submerged  body  in  motion  are  developed  in  relation  to  the  shape  of  the  vehicle  and  are 
more  conveniently  expressed  in  relation  to  the  body  axes.  Equating  the  applied  force 
vector  to  the  acceleration,  results  in, 

F  =  m/v  +  wxpG  +  wxwxpG+coxpGj  (213) 

The  applied  force  vector  is  composed  of  gravitational  (weight)  and  hydrostatic  (buoyancy) 
forces  together  with  any  hydrodynamic  forces  arising  from  relative  motion  between  the 
vehicle  and  the  ocean  water  particles  These  will  be  discussed  in  further  detail  in  following 
sections 

7.  Rotational  Equations  Of  Motion 

To  develop  the  rotational  equations  of  motion,  the  sum  of  applied  moments 
about  the  vehicle's  center  of  mass  is  equated  to  the  rate  of  change  of  angular  momentum  of 
the  vehicle  about  its  center  of  mass.  This  equating  will  provide  the  necessary  equations  of 
motion  for  the  remaining  three  degrees  of  freedom  In  the  practical  case  of  marine 
vehicles,  however,  the  statement  just  made  is  modified  slightly  because  it  is  much  more 
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difficult  to  assess  the  vehicle's  mass  moments  of  inertia  about  its  center  of  gravity  (CG) 
As  the  CG  changes  with  loading;  it  becomes  simpler  to  evaluate  the  mass  moments  of 
inertia  about  the  body  fixed  frame  which  lies  along  the  axes  of  symmetry  of  the  vehicle  in 
most  cases.  The  mass  moments  of  inertia  for  the  vehicle  to  be  computed  are; 


1  = 


xx  xy  x 

y*        >y        y 


i. 


i„ 


(214) 


The  angular  momentum  of  the  body  is  thus; 

#0  =  /0w,  (2.15) 

resulting  in  the  total  applied  moments  about  the  origin  given  by, 

M0=4+PgxMg/  (2.16) 

Since  the  rate  of  change  of  angular  momentum  is  given  by, 

fl0  =  I06)  +  (bxH0,  (2.17) 

and  the  acceleration  of  the  global  position  vector  is  given  by, 

^o  =  v-  +  coxv,  (2  18) 

then  the  rotational  equation  of  motion  in  vector  terms  is  given  by, 

M0  =I0(b  +  (bx(I0(&)  +  m{pGxv  +  pc  xibxv}.  (2  19) 

The  applied  moments  about  the  body  fixed  frame  arise  from  a  static  balance  of 
weight  and  buoyancy  effects  to  achieve  the  proper  trim  and  heel,  and  hydrodynamic 
moments  from  the  forces  applied  through  appendages  such  as  control  fins,  hydrodynamic 
effects  from  waves,  and  hydrodynamic  effects  from  relative  motion  between  the  vehicle 
and  the  water. 

At  this  point,  there  are  three  translational  equations  obtained  from  equation 
2.13,  three  rotational  equations  obtained  from  equation  2.19,  and  six  unknown  velocities 
(wand  v).  In  itself  this  is  aconsistent  set  of  dynamic  equations  if  the  weight  and  buoyancy 
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terms  are  always  self  canceling.  However,  with  weight  and  buoyancy  being  applied  in  a 
global  vertical  direction,  and  also  being  applied  at  different  locations,  they  represent  forces 
that  are  dependent  on  the  attitude  of  the  vehicle.  Notice  that  without  weight  and  buoyancy 
forces,  the  six  consistent  equations  expressed  in  the  body  fixed  frame  of  the  vehicle  could 
be  solved  for  the  vehicle's  velocity  (given  applied  force  descriptions).  With  weight  and 
buoyancy  acting,  we  need  additional  equations  (constraints)  that  will  link  vehicle  attitude 
to  motion  so  that  the  combined  set  of  equations  can  be  solvable.  The  constraints  to  be  used 
are  developed  using  the  Euler  angles  which  are  a  general  set  of  angles  used  to  define  the 
attitude  of  the  rigid  body.  From  these  angles  a  relationship  between  rate  of  change  of 
attitude  and  the  body  rotational  rates  given  earlier  as  the  angular  velocity  vector 
63  =  / p.q.r  J ',  can  be  developed 

8.      Incorporation  Of  Vertical  Forces  Into  The  Equation  Of  Motion 

The  weight  and  buoyant  forces  that  act  at  the  centers  of  gravity  and  buoyancy 
must  be  defined  from  static  analyses.  For  submerged  bodies  the  weight  and  buoyancy  force 
vectors  do  not  change  with  vehicle  attitude  For  a  surface  ship,  B  will  change  with 
attitude,  and  the  righting  moments  must  be  computed  from  Naval  Architectural 
considerations  Assuming  that  weight  and  buoyancy  are  fixed  in  relation  to  the  body  fixed 
frame,  W  =  0/  +  0 J  +  (mg)K,  and  B  =  0l  +  0J -(pV )K .  Since  the  weight  and  buoyancy 
terms  in  the  applied  forces  act  in  the  global  vertical  direction,  they  must  be  transformed 
into  components  in  the  vehicle  fixed  frame  before  they  can  be  added  into  the  equations  of 
motion  Returning  to  equation  2  4,  we  therefore  find  the  components  along  the  vehicle 
fixed  frame  to  be  the  third  column  of  that  transformation  matrix  The  net  vertical  force 
components  become. 
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l=(W-B) 


-sinQ 
cos  Q  sin  § 
cos 0  cos  (f) 


(2.20) 


The  weight  portion  of  the  vertical  force  acts  at  the  center  of  gravity  of  the  vehicle  The 
buoyancy  portion  of  the  vertical  force  acts  at  the  center  of  buoyancy.  Because  these  forces 
act  at  positions  away  from  the  body  center  they  exhibit  a  moment  about  the  body  center 
given  by, 


m=WpGx 


-sin® 

-  sin  8 

cosQsinty 

-BpBx 

COS  05/>7(j) 

cosQcosty 

cos  Q  cos  § 

(2.21) 


This  moment  will  not  be  zero  even  if  W  and  B  are  identical  because  it  is  not  usually  the 
case  that  pG  and  pg  are  collocated.  For  static  stability  it  is  advisable  to  locate  the  center  of 
gravity  below  the  center  of  buoyancy  The  total  vertical  force  vector  can  be  written  as, 


Fc  = 


fa 


Vflr 


(2.22) 


and  is  added  negatively  to  the  left  hand  side  of  the  equations  of  motion 

9.      Development  Of  The  Full  Six  Degrees  Of  Freedom  Non-Linear  Equations 

Of  Motion  For  A  Marine  Vehicle 

Define  a  vector  x  of  vehicle  body  frame  velocities  to  be,  x  =  fu,v,w,p,q,rj , 

and  a  vector  z  of  global  positions  to  be  z  =  [X,Y,Z,$,Q,\\f  ]\  then  considering  M  as  a 

6x6  mass  matrix  including  translational  and  rotational  inertial  elements,  the  equations  of 

motion  can  be  written  in  the  following  vector  form, 

MZ  +  f(x)+Fg(z)  =  Fh  (2.23) 

and, 

i  +  g(x,z)  =  0.  (2.24) 
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Therefore  with  suitable  knowledge  of  the  excitation  force  and  moment  loads, 
Fh,  as  a  function  of  time  and/or  vessel  motion,  a  solution  for  the  vehicle's  dynamics  can  be 
obtained  A  more  detailed  insight  into  the  development  of  the  twelve  differential  equations, 
in  first  order  form  given  by  the  foregoing  analysis  shows, 


m{v  +  CO  x  pG }  +  m((b  x(bxpG+(bxv)  +  fG(z)  = 


X. 


and. 


(2.25) 


I0U)  +  mfpG  xv}  +  (&x(I0(b)  +  mfpG  xaixv}  +  mg(z)  = 

It  helps  here  to  define  the  cross  product  coefficient  matrix  so  that, 

co  x  pG  =  -pG  x  co  =  -/"cros(pG  )7co 
where, 


(2.26) 


(2.27) 


fCTOS(pG)J  = 


0     zo    -yG 

-zG       0        xG 


(2.28) 


yG    ~xg     ° 

Now  collecting  the  mass  matrix  coefficients  into  a  6x6  matrix  including  the  inertia  cross 
coupling  effects,  results  in  the  following  equation, 
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M  = 


m 

0     0 

0 

~G           "^ 

0 

m     0 

/w 

~ZG 

0        xG 

r 

0 

0 

0     m 

_  f 

*G 

•VC   " 

yG 

-xG       0 

7 
**G 

0      -xG 

/„ 

vy           vz 

.- 

?G 

*G 

0 

/« 

K  I*_ 

m    zr 


(2.29) 


The  remaining  terms  on  the  left  hand  side  of  the  equations  of  motion  arising  from  the 
centripetal  and  coriolis  accelerations  become, 


f(x)  = 


m((£>  x  cb  x  pG  +  co  x  v) 
co x (I0w)  +  mfpG  x(bxvj 


(2.30) 


The  double  vector  cross  products  are  nonlinear  in  the  primary  velocity  variables  and  hence 
the  need  for  the  nonlinear  functional,  /(•)  The  reader  can  perform  the  indicated 
manipulations  to  express  individual  equations  within  the  set  if  so  desired. 

The  components  of  the  hydrodynamic  and  external  forces  and  moments  acting 
on  the  vehicle  body  are  separated  in  the  above  analysis  into  six  components  each  acting 
along  the  vehicle  body  fixed  coordinate  axes  and  form  the  total  vector  of  forces  and 

moments  as, 

Fh(t)  =  [Xf(t),Yf(t),Zf(t),Kf(t)Mf(t),Nf(t)]\  (2.31) 

where  the  vector  components  in  order  refer  to  the  surge1  sway,  heave  forces,  and  the  roll, 
pitch,  and  yaw  moments  respectively. 

The  long  form  of  the  six  degrees  of  freedom  equations  of  motion  can  be  written 
as  follows, 


m[u-vr  +  wq-xG(q2+r)  +  yG(pq-f)  +  zG(pr  +  q)]  = 
-(W-B)sin9+Xf 


(2.32) 
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m[v  +  ur-wp  +  xG(pq  +  r)-yG(p2+r)  +  zG(qr-p)]: 
(W-B)cos§sm§  +  Yf 

m[\i-uq  +  vp  +  xG(pr-q)  +  yG(qr  +  p)-zG(p2+q2)] 
(W-B)cosQcosQ>  +  Zf 


(2.33) 


(234) 


IJ  +  (I:-Iy)qr  +  IJpr-q)-IJq2-r2)-IJpq  +  r)  + 

m[yG(™  ~uci  +  VP)  -zG(v  +  ur-  wp)J  =  (yGW-yBB)cosQcos<^  (2.35) 

-(zG  W  -  zG  B)  cos  6  sin  <f>  +  Kf 

IJ  +  (Ix-IJpr-L(<ir  +  p)  +  IJpq-f)  +  IJp2-r2)- 

m[xG(w  -uq  +  vp)-zG(u  -vr  +  wq)]  =  -(xGW  -  xBB)cosQcos<\>  (2.36) 

-(zGW-zBB)smQ  +  Mf 

I:>-  +  (ly-Ix)Pq-IJp2-<l2)-IyJpr  +  q)  +  lJqr-p)  + 

m[xG(v  +  ur-wp)-yG(u-vr  +  wp)J  =  (xGW-xBB)cosQsin§  (2.37) 

+(yGW-yBB)sind+Nf 

while  the  kinematic  relations  for  the  vehicle  rate  of  change  of  attitude  and  motion  over  the 
ocean  bottom  require  equations  2  1 1  and  2.8. 

10.    Adaptation  Of  The  Non-Linear  Equations  Of  Motion  To  The  Vertical 

Plane 

Restricting  the  motions  of  the  vehicle  to  the  vertical  (dive)  plane,  the  only 
significant  motions  that  must  be  incorporated  to  effectively  model  the  vehicle  in  the  dive 
plane  are,  the  surge  velocity  (w),  the  heave  velocity  (w),  the  pitch  velocity  (q),  the  pitch 
angle  (6)  and  the  global  depth  position  (Z)  This  restriction  simplifies  the  twelve 
previously  developed  equations  to  a  system  of  four  non-linear  equations  of  motion,  which 
are; 

m(w-uq-zGq2-xGq)  =  Zf  (2.38) 


l: 


Iyq  +  mzGwq-mxG(w-uq)  =  Mf  (2.39) 

0  =  4  (2  40) 

Z  =  -U  57/70  +  wcosQ  (2  41) 


where, 


Zf  =  Zqq  +  Z.w  +  Zquq  +  Zwuw  - 

1  nose  t     —      \^  (^  42) 
^p  J  Q/W*'         ,  dx  +  (W-B)cosd+u\Zs8s  +  ZSb8b) 

2  L  lH'"^l 


and, 


A/,= 


,.  =  M .a  +  M  .w+ M  uq  + M  uw 
f         q*        w  q  *        w 


1    7Ci(x)(H'    *^)  xdx 
-±p\  \w-xq\  ,        (2.43) 

*"      tail 

-(x„W - x DB)cosd-(z„W - zDB)s\nd+u2(Ms  8S  +  M5  8b) 
results  from  expanding  the  Zf  and  the  M7  terms  from  (2.3 1)  in  a  first  order  Taylor  series 

expansion  and  incorporating  both  hydrostatic  and  fluid  drag  forces. 

Equations  2.38,  2.39,  2.40  and  2.41  can  be  linearized  for  a  level  flight  path  when 
the  dive  plane  angle  is  zero,  i.e.  5o=0.  By  setting  all  time  derivatives  to  zero,  and 
neglecting  for  the  moment  the  hydrodynamic  drag  terms  the  following  are  obtained, 

Zjuw  +  (W-B)cos$  =  0  (2  44) 

Mjiw-(xcW-xBB)cosQ-(zGW-zBB)sine  =  0  (2  45) 

q  =  0  (2  46) 

-usinQ  +  wcosQ  =  0.  (2  47) 
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If  the  assumption  is  made  that  the  vehicle  is  neutrally  buoyant,  then  it  can  be  said  that 
xG  =  xB  and  W  -B  From  this  equations  2  44  and  2.45  can  be  reduced  to, 

Zjuw  =  0  (2.48) 

Mwuw-(zG-zB)BsinQ  =  0,  (2.49) 

which  yield  the  nominal  position  wo  =qo=0  and  sinQ0  =  0,  resulting  in  the  0O  solution  as 
either  90  =  0  or  0O  =  n   Choosing  to  linearize  around  the  nominal  point  0O  =  0  results  in; 

q2  =  2qoq  =  0 
wq=woq  +  qow  =  0 
sinQ  =  cosQoe  =  Q 
wcosQ  =  (-MosmQo)Q  +  (cosQo)w  =  w . 

The  linear  equations  of  motion  are  then  written  as; 

(m  -  Zw  )w  -  (Zq  +  mxG  )q  =  Zwuw  +  (Zq  +  m)uq  +  Z5w25 

(Iy  -Mq)q- (MH  +  mxG )w  =  Mwuw  +  (Mq  - mxG )uq - (zG  -zB)m  +  M5w:8 

6  =  <7 

Z  =  -uQ  +  w 
As  equations  2  42  and  2.43  show,  both  Z§  and  Mg  are  a  linear  combination  of 

the  respective  stern  and  bow  hydrodynamic  control  surface  coefficients  and  the  respective 

input  value  of  5    This  makes  the  system  of  equations  as  a  multiple  input  system.  To 

reduce  this  system  into  a  single  input  system  the  linear  combination  of  control  inputs  will 

be  modified  into  the  following  form, 

Zh=(Zb  +aZ6  ;.  (2.58) 

s  s 

This  will  allow  a  single  input  5  to  control  both  stern  planes  and  bow  planes,  and  will 
cause  the  bow  planes  to  be  slaved  to  the  stern  planes  This  technique  is  known  as  dual 
control  The  value  (X  will  range  from  -1  to  1  The  selection  of  the  value  of  (X  will  allow 
the  planes  to  operate  as  desired  for  the  particular  maneuvering  condition,  i.e.,  a  =  0  for  no 
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(2.50) 

(2.51) 

(2.52) 

(2.53) 

(2.54) 

(2.55) 

(2.56) 

(2.57) 

bow  plane  control,  a  =  -1  for  bow  plane  and  stern  plane  control  opposed  to  each  other, 
yielding  the  maximum  pitch  moment,  and  a  =  1  for  bow  and  stern  plane  control  in  the 
same  direction,  yielding  the  maximum  heave  force 

In  state  space  form  these  equations  can  be  written  as. 


"e" 

0 

0 

1 

0" 

"e" 

0 

<7 

= 

a\  yzGB 

a2lZGB 

a2]u 

a22u 

0 
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+ 

bxir 
b2u2 

Z 

-u 

1 

0 

0 

z 

0 

6, 


where  the  coefficients  at]  and  bt  are  given  by; 

Dv=(m-ZJ(Iy-Mq)-(mxG  +  ZqXmxG+MJ 

auDv  =  (Iy-Mq)Zw  +  (mx0+Zq)Mw 
avDv=(Iy-Mq)(m+Zq)  +  (mxG  +  Zq)(Mq-mxG) 
a]}Dv  =  -(mxG  +  Zq)W 
b,Dv  =  (Iy  -  Mq)Zb  +  (.mxG  +  Zq)Mh 

a2]Dv  =  (m-  ZJMw  +  (mxG  +  MJZw 
a22Dv  =  (m-ZJ(Mq  -mxG)  +  (mxG  +  MJ(m  +  Zq) 

a2iDv  =  -(m-ZJW 
b2Dv  =  (m-  Zw)Mh  +  (mxG  +  Mw)Z6 
and  zGB  =  zG-zB  is  the  metacentric  height. 

This  linear  system  of  equations  becomes  the  basis  for  the  control 
the  following  section 


(2.59) 


(2.60) 

(261) 

(262) 

(263) 

(264) 

(265) 

(2.66) 

(267) 

(268) 

design  in 

B.     CONTROL  LAW  DESIGN 

1.      Introduction 

The  control  design  problem  can  be  stated  as  follows:  Given  the  system, 
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x  =  Ax  +  Bb, 


(269) 


where  the  state  vector  equation  is, 


x  = 


(2.70) 


how  do  we  find  5,  such  that  the  system  will  behave  as  desired.  The  type  of  control  that  is 
of  interest  in  this  problem  is  closed  loop  control,  where  5  is  a  function  of  the  state  x 
Since  the  state  x  is  used  to  determine  the  control  effort  h(x)  it  is  called  feedback  control. 

2.      Pole  Placement 

A  linear  full  state  feedback  control  law  is  introduced  in  the  form 

b  =  -Kx,  (2.71) 

where  K  is  the  feedback  gain  vector  to  be  determined  such  that  the  closed  loop  system  of 
equations  2.69  and  2  71  has  the  desired  system  dynamics.  Substituting  equation  2.71  into 
equation  2  69  yields, 

x  =  (A-BK)x.  (2.72) 

The  actual  characteristic  equation  of  this  closed  loop  system  is  given  by 

det{A-BK-sl]  =  0.  (2.73) 

The  gain  vector  K  can  be  chosen  such  that  the  actual  characteristic  equation  assumes  any 
desired  set  of  eigenvalues.  If  the  desired  locations  of  the  closed  loop  poles  are  chosen  at 
s  =  st  for  /-l, ...,«,  the  desired  characteristic  equation  becomes 

(s-sj(s-s:)...(s-sj  =  0  (2.74) 

The  required  values  of  K  are  obtained  by  matching  coefficients  in  the  two  polynomials  of 
the  actual  and  desired  characteristic  equations 
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Consider    the    previously    developed    linear    state    matrix    equation    2.59 
Substituting  equation  2.71  into  equation  2.59  results  in  the  closed  loop  system 


a 


-b]u2k]     auu-b^u'k2     a^u-b^urk^     -b]u~k4 
azizGB  ~  ^yK     az\u  ~  bzu2k2     a22u  -  b2u2ki     -b2u2k4 


13*GB 


-u 


1  0  0 

The  characteristic  equation  of  the  closed  loop  system  is  given  by 


(2.75) 


det 


-s  0  1 

anzGB  ~ b\u K  auu-b{n2k2-s  anu-bxu2k^ 

a2JzGB  -  b2u~k\  a2]u-b2u2k2  a22u  -  b2u'k3  -  s 

-u  1  0 


0 
-b^u2k4 
-b2u2k4 

-s 


=  0, 


which  after  some  algebra  reduces  to 


where , 


s4  +  (  A2k2  +  A,k}  -  £,  >3  +  (-B,k,  -  B2k2  -  B3k,  -  B4k4  -  E2  Js2 
+  (-C]k]-C2k2-C4k4-E,)s+-(-D]k4-D2kJ  =  0 

A,  =  -B4  =  b}u~ 
Aj  =  -B]  =  b2u2 


B2  =(a22bx  -aub2)u3 

C2  =  D,=(a2&-anb2)zGBu2 


C4  -  (a22bx  +b2  -avb2)iil 
D2  =(aub2  -a^bju4 
Ex  =(an  +a22)u 

E2  ~  a23ZGB  +  fai:a:i  _aiia22>2 


(2.76) 


(2.77) 

(2.78) 

(2.79) 
(2.80) 
(2.81) 
(2.82) 
(2.83) 
(284) 
(285) 
(2.86) 
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£3  =  (ana2]  -aua2iJzGBu  (2.87) 

Now  if  the  closed  loop  poles  are  placed  at  -/?,,  -p2,  -  p3,  and  -p4,  then  the 
desired  characteristic  equation  is, 

(s  +  P])(s  +  p:Xs  +  p3Xs  +  pJ  =  0,  (2.88) 


or 


with. 


54  +  a]53  +  a:5:+a35  +  a4  =  0  (2  89) 


a,  =A+A  +  A+A  (2.90) 

<*:  =AA+AA+AA+AA+AA+AA  (291) 

a3  =  AAA  +  AAA  +  AAA  +  AAA  (2.92) 

«4=AAAA  (293) 

The  control  gains  can  now  be  computed  by  equating  coefficients  of  the  actual  and  desired 

characteristic  equations, 

A2k2  +  Aik1  =  -a]-E]  (2.94) 

#,*,  +  B2k2  +  B,k,  +  B4Jc4  =a2  +  E2  (2.95) 

C,*,  +  C2k2  +  C4k4  =  a3  +  £3  (2.96) 

(I\+Ch)kA=a.A.  (2.97) 

Now  that  a  method  for  computing  the  gains  of  a  controllable  single  input  system 

that  will  place  the  poles  at  any  desired  location  has  been  developed,  the  selection  of  pole 

location  must  be  accomplished. 
3.      Pole  Location  Selection 

In  a  typical  second  order  system  control  law  design  the  transient  response 

specifications  will  be  given  This  results  in  an  allowable  region  in  the  s-plane  where  the 

desired  location  of  the  poles  can  be  obtained.  For  higher  order  systems  the  concept  of 

dominant  roots  can  be  employed.  In  selecting  poles  for  a  physical  system  it  is  necessary  to 
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look  at  the  physics  of  the  system.  If  the  poles  are  specified  too  negative,  a  very  small  time 
constant  for  the  control  system  will  result,  and  the  physical  system  may  not  be  able  to 
react  that  fast  The  control  law  u  =  -Kx,  implies  that  for  a  given  state  x  the  larger  the 
gain,  the  larger  the  required  control  input  In  practice  there  are  limits  typically  placed  on  u 
(actuator  size  and  saturation).  Occasional  control  saturation  is  not  serious  and  may  even 
be  desired  A  system  that  never  saturates  is  in  all  likelihood  over  designed 

Considering  the  control  law  design  to  stabilize  the  submarine  to  a  level  flight 
path  at  0  =  0  it  will  be  required  that  the  submarine  return  to  level  flight,  after  a  small 
disturbance  in  6  or  Z,  within  the  time  it  takes  for  the  vehicle  to  travel  three  ship  lengths 
Since  the  submarine  is  about  14  feet  long  and  it  travels  at  5  ft/sec,  the  required  recovery 
time  is  about  10  seconds.  This  means  that  the  time  constant  is  about  3  seconds,  and  the 
closed  loop  poles  should  be  placed  at  approximately  -0.3. 

Control  design  using  pole  placement  is  very  easy  using  MATLAB  The 
appropriate  MATLAB  command  is  place,  which  accepts  as  inputs  the  A  and  B  matrices 
along  with  a  vector  of  the  desired  closed  loop  poles,  and  returns  the  gain  vector  K.  Using 
equation  2.59  with  a  nominal  speed  of  5  ft/sec  and  the  vector 
p  =  f-0.3  -0.31  -0.32  - 0.337,  (place  does  not  like  poles  in  the  exact  same  location), 
the  gain  vector  K  was  calculated  using  MATLAB  and  by  simultaneously  solving  equations 
2.94  through  2.97  yielding  the  same  results  Substituting  the  gain  vector  and  state  variable 
vector  into  equation  2.71  results  in  a  control  law  of 

5  =  -f-O.9917e-0.8333w-O.6O26<7  +  0.O351Z/  (2  98) 

Using  this  control  law  along  with  equations  2.38  through  2  43  a  simulation 
program  was  developed  to  investigate  the  vehicle's  response  to  initial  disturbances,  the 
results.of  which  will  be  presented  in  the  following  chapter 
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m.  PROBLEM  IDENTIFICATION 


A.     VEHICLE  SIMULATIONS 

With  the  control  law  developed  and  the  vertical  plane  equations  of  motion  defined,  a 
vehicle  simulation  program  was  developed  (see  Appendix  A).  This  program  was  used  to 
simulate  the  vehicle's  response  at  a  variety  of  speeds  and  initial  disturbances.  The 
simulations  were  conducted  with  several  simplifications  applied  to  the  vertical  plane 
equations  of  motion  (E.O.M.),  equations  2.38  through  2.43.  The  simplifications  or 
constraints  applied  to  the  E.O.M.  were  , 

1 )  to  consider  the  body  drag  forces  negligible, 

2)  to  consider  the  vehicle  to  be  neutrally  buoyant,  i.e.,  (W=B), 

3)  to  assume  that  the  locations  of  the  center  of  gravity  and  center  of  buoyancy, 
in  the  X-Y  plane,  are  coincidental,  and 

4)  that  the  rudder  is  restricted  to  ±  23  degrees. 

These  assumption  resulted  in  the  reduction  of  equations  2.38  through  2.43   to  the 
following  system  of  equations, 


1 

0 

0 

0 

e 

0 

(m-ZJ 

-^ 

0 

vi' 

0 

-Mw 

h  -  Mq 

0 

<7 

0 

0 

0 

1 

Z 

m(uq  +  zGq2  )  +  Zquq  +  Zwuw  +  u% 
m(-zGwq)+  M  uq  +  Mwiw  -  BzGB  sinQ  +  uzMb 


-usinQ  +  wcosQ 


(3.1) 


forming  the  basis  of  the  initial  vehicle  simulations 
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The  vehicle's  response  to  small  disturbances,  as  stated  earlier,  was  investigated  over  a 
wide  range  of  speeds.  A  small  sample  of  these  simulations  is  shown  in  Figures  3-1  through 
3-3. 

As  Figure  3-1  depicts,  the  vehicle  returns  to  level  flight  in  approximately  10  seconds. 
This  response  is  as  expected  based  on  the  control  system  design  of  Chapter  II.  The 
simulation  results  represented  in  Figure  3-2  also  show  the  vehicle  steadying  out  in  a  level 
flight  at  the  desired  depth,  however  the  vehicle  requires  a  significant  amount  of  time  for 
this  to  occur.  The  reason  for  this  extended  amount  of  time  is  that  the  control  law  gains 
were  set  for  the  nominal  operating  speed  of  5  fps,  and  the  vehicle  is  being  simulated  at 
one-third  of  that  speed  This  results  in  a  very  slow  return  to  ordered  position  caused  by 
the  reduced  hydrodynamic  effects  at  this  speed 

The  final  plot  of  the  sample  vehicle  simulations,  Figure  3-3,  illustrates  a  problem  with 
the  vehicle  stability.  As  can  be  seen,  the  vehicle  does  not  return  to  the  desired  position,  but 
instead  steadies  out  at  a  pitch  angle  of  eight  degrees,  with  a  steady  state  dive  plane  value 
of  23  degrees  This  same  type  of  response,  but  with  different  pitch  angle  values,  was 
observed  in  all  the  simulations  conducted  below  approximately  19  fps  These  results 
indicate  that  there  is  a  critical  point,  or  speed  at  which  the  vehicle  stability  is  affected 
Determining  this  point  will  be  the  basis  for  the  remainder  of  this  chapter. 

B.     CRITICAL  SPEED  IDENTIFICATION 
1.      Eigenvalue  Analysis 

Consider  the  nonlinear  system  of  state  equations  3.1,  which  can  be  written  in  the 
form, 

x  =  f(x).  (32) 

It  is  known  that  the  equilibrium  points,  xo,  of  the  system  are  defined  by 
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THET A/DELTA  vs.  TIME 
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Figure  3-1.  Vehicle  Response  For  a  Nominal  Operating  Speed  of  5  fps. 
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THET A/DELTA  vs.  TIME 
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Figure  3-2.  Vehicle  Response  at  2.0  fps. 
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THET A/DELTA  vs.  TIME 
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f(xj  =  0.  (3.3) 

This  is  a  nonlinear  system  of  algebraic  equations  and  it  may  have  multiple  solutions  in  x0, 
which  means  that  the  nonlinear  system  of  equations  may  have  more  than  one  position  of 
static  equilibrium.  If  one  equilibrium  value  x0,  is  selected,  the  stability  properties  can  be 
established  by  linearization.  The  linearization  converts  equation  3.2  into  the  form 

x  =  Ax,  (3.4) 

where  A  is  the  Jacobian  matrix  offfx)  evaluated  at  the  equilibrium  point  x0, 

a/ 


A  = 
dx 


(3.5) 


and  the  state  x  has  been  redefined  to  designate  small  deviations  from  the  equilibrium  x0, 

x^x-x0  (3.6) 

As  long  as  all  eigenvalues  of  A  have  negative  real  parts,  the  linear  system  will  be  stable 
This  means  that  the  equilibrium  x0  will  be  stable  for  the  nonlinear  system  as  well.  This 
statement  is  nothing  more  than  Lyapunov's  linearization  technique 

The  question  that  must  be  answered  is,  what  happens  if  one  real  eigenvalue  of 
the  linearized  A  matrix  is  zero7  The  interesting  case  here  is  when  the  rest  of  the 
eigenvalues  have  all  negative  real  parts,  otherwise  x0  is  unstable  and  the  problem  is  solved. 
If  the  case  of  a  zero  eigenvalue  appears  too  specialized  to  be  of  any  practical  use,  consider 
this:  Assume  that/fx)  depends  on  one  physical  parameter  and  that  physical  parameter  is 
allowed  to  vary  over  some  range.  Then  clearly  A  will  depend  on  that  parameter  and  as  the 
parameter  varies,  it  is  possible  that  one  real  eigenvalue  of  A  will  become  zero  for  a  specific 
value  of  the  parameter.  The  problem  is  then  to  establish  the  dynamics  of  the  nonlinear 
system  as  one  real  eigenvalue  of  A  crosses  zero,  i.e.,  goes  from  negative  to  positive  As 
the  solution  evolves  in  time,  things  are  interesting  only  along  the  direction  of  the 
eigenvector  that  corresponds  to  the  critical  eigenvalue.  Along  the  rest  of  the  directions  in 
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the  state  space,  everything  should  converge  back  to  the  equilibrium,  remember  it  was 
assumed  that  all  remaining  eigenvalues  of  A  have  real  negative  parts.  Although,  strictly 
speaking,  this  is  a  true  statement  for  linear  systems,  there  are  technical  reasons  that  force  it 
to  be  true  for  nonlinear  systems  as  well,  the  only  difference  is  that  the  corresponding 
directions  in  the  state  space  are  curved  instead  of  straight. 

By  taking  the  previously  developed  characteristic  equation  2.77,  and 
remembering  that  the  roots  of  this  equation  are  the  eigenvalues  of  interest,  it  is  easy  to  see 
that  if  a4  of  equation  2.93,  is  set  equal  to  zero,  then  one  eigenvalue  will  be  zero  If  this  is 
done,  and  recalling  that  k4  holds  a  non-zero  value  equation  2.97  reduces  to 

(°\  A  "  a:  A  V  +  zgb  (a\  A  ~  azA  )  =  °  (3  7) 

Since  all  parameters  of  equation  3.7  have  a  fixed  value  with  the  exception  of  u,  there  must 
be  some  value  of  u  that  causes  the  linear  A  matrix  to  become  unstable.  Recognizing  this 
fact  will  allow  equation  3  7  to  be  solved  in  terms  of  u  yielding 


//  = 


ZGS(a?}>\-a\A) 


-An 


(3.8) 


(aub2-a2fo) 
Using  equations  2  60  through  2  68,  equation  3  8  can  be  simplified  into  the  following  form, 


u  = 


"    7  R 


-il/2 


(MWZ8-ZWM8 

Substituting  the  appropriate  values  from  Appendix  B,  equation  3.9  can  be  solved  for  the 
value  of  u  that  causes  the  systems  to  become  unstable  This  results  in  a  value  of  u  = 
1  8979  fps  using  an  a  =  0.  This  value  corresponds  very  closely  to  the  value  of  1  9  fps 
obtained  by  vehicle  simulations. 


(3.9) 
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As  a  check  to  this  solution,  a  MATLAB  program  was  written  to  compute  the 
closed  loop  eigenvalues  of  the  matrix  equation  2.59  as  a  function  of  the  speed  u 
Figure  3-4  shows  a  plot  of  the  real  parts  of  the  eigenvalues  versus  speed  It  can  be  seen 
that  the  real  part  of  the  eigenvalue,  represented  by  the  dash-dot  line,  becomes  positive  at 
speeds  less  than  approximately  19  fps  supporting  the  earlier  findings.  (For  clarity 
purposes,  the  values  of  the  eigenvalue  represented  by  the  dash-dot  line  were  scaled  by  a 
factor  often.) 

2.      Steady  State  Analysis 

To  determine  the  existence  of  multiple  steady  state  solutions,  for  the  system 
represented  by  equation  2.59,  it  is  necessary  to  perform  a  steady  state  analysis  on  the 
system  of  equations  which  models  the  vehicle  dynamics  Using  matrix  equation  3.1,  and 
setting  all  time  derivatives  equal  to  zero,  results  in  the  following  set  of  equations; 

q  =  0,  (3  10) 

Zquq  +  Zwuw  +  u2Z&b  =  m(-uq-zGq2),  (3  11) 

Mquq  +  Mwuw  -  BzGB  sinQ  +  u~M&5  =  mzGwq  ,  (3.12) 

-usind  +  wcosQ  =  0  (3.13) 

Simplifying  these  four  equations  yields  the  following  two  equations  in  terms  of  6 ,  u,  and 
the  physical  parameters  of  the  vehicle, 

ZytanQ  +  u2Zfi  =  0  (3  14) 

Mwu2  tanQ  -  BzGB  smQ  +  u2Msb  =  0 .  (3  15) 

Now  if  equation  3.14  is  multiplied  by  Ms,  equation  3.15  is  multiplied  by  Z5  and  the  two 

resulting  equations  are  set  equal  to  each  other  and  simplified,  the  following  single  equation 
is  obtained, 

[(ZwMb  - MJb)u2  +  ZbBzGB cosQJsinQ  =  0.  (3.16) 
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Real  Part  of  Eigenvalues  vs.  Speed 
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Figure  3-4.  Real  Part  of  Closed  Loop  System  Eigenvalues  as  a  Function  of  Speed. 
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Through  inspection  it  can  be  seen  that  equation  3  16  could  have  multiple 
solutions  in  6,  besides  the  trivial  solution  of  6  =  0.  So,  if  this  equation  is  rearranged  and 
solved  for  9*0,  solutions  occur  when 

hVkz^K,  (317) 

7  R- 

Therefore,  the  steady  state  value  of  6  is  represented  by 

u2(MwZ6-ZwMj 


6„  =  cos' 


ZbBzGB 


(3.18) 


which  may  have  multiple  solutions  based  on  the  value  of  the  speed  u.  As  discussed 
previously,  these  multiple  solutions  were  evident  in  the  simulations  conducted  earlier  in 
this  chapter 

Knowing  that  the  maximum  value  for  cosQ  is  equal  to  one,  the  right-hand  side 
of  equation  3  17  will  be  true  for  values  less  than  or  equal  to  one.  If  this  constraint  is 
imposed  on  3  17,  then  the  equation  may  be  manipulated  into  the  following  form 

uz< 5^ .  (3  19) 

(MJLh-ZMO 

Rearranging  this  equation  and  solving  for  the  upper  limiting  case  yields  the  same  results  as 
equation  3.9.  Again  this  supports  the  critical  speed  values  discussed  earlier 

Equation  3.19  can  be  expressed  in  a  non-dimensional  Froude  number  form  by 
dividing  the  equation  by  both  the  gravitational  constant  g  and  the  physical  parameter  zGB 
and  taking  the  square  root  of  the  result.  This  will  convert  equation  3.19  into  the  non- 
dimensional  form  given  below, 


Fn=l^-<    5* (3.20) 
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With  the  steady  state  9  solution  derived,  equations  for  both  the  steady  state 
value  of  5  and  steady  state  value  of  Z  can  be  determined  Using  equation  3  14  and 
substituting  the  steady  state  value  of  6   will  allow  the  steady  state  5   equation  to  be 

obtained. 

o„  =  ^to»6„  (3.21) 

To  achieve  the  steady  state  solution  of  Z,  begin  by  writing  equation  2.71  in 
expanded  general  form  as 

^  =  -(k]e  +  k2w  +  k3q  +  k4Z).  (3.22) 

Then  by  applying  the  same  conditions  to  equation  3.22  that  resulted  in  equations  3  14  and 
3  1 5,  and  using  the  previously  developed  steady  state  equations  for  0  and  5 ,  an  equation 
of  the  form 

-7 

— s-tanQ.  +  kfl.  +  kptanQ. 

Z„=^ (323) 

-K 

is  obtained. 

The  loss  of  stability  of  an  equilibrium  and  the  generation  of  additional 
equilibrium  states,  is  called  a  pitchfork  bifurcation  and  is  very  common  in  nature  The 
buckling  of  a  beam  is  one  such  example.  Using  equations  3.18,  3  21,  and  3.23  as  a  basis 
for  another  MA  TLAB  program,  the  steady  state  solutions  of  0 ,  5 ,  and  Z  versus  Froude 
number  were  investigated  with  the  results  displayed  in  Figures  3-5  through  3-7.  These 
three  plots  are  referred  to  as  supercritical  pitchforks,  so  named  because  upon  the  loss  of 
stability  of  the  trivial  equilibrium  the  additional  nearby  equilibrium  states  are  stable 
Graphically,  this  can  be  represented  in  Figures  3-5  through  3-7,  where  the  solid  curves 
represent  stable  and  dotted  curves  represent  unstable  equilibria  [Ref  17].  These  three 
figures  also  demonstrate  the  control  input  saturation.  Upon  control  surface  saturation,  the 

36 


steady  state  equations  for  0 ,  5  and  Z  are  no  longer  valid  due  to  the  limit  placed  on  the 
maximum  plane  angle  The  steady  state  equations  that  hold  for  this  region  of  saturation 
are  given  below, 

5„=±  22.9°  (0.4  radians),  (3.24) 

"    =  ^,  (3.25) 


and, 


6„  =  sin 


ZGBB 


(326) 


As  can  be  seen  in  Figures  3-5  and  3-6,  at  an  approximate  Froude  number  of  1.05  the 
control  surfaces  saturate  resulting  in  the  linear  solutions  displayed  below  that  saturation 
Froude  number  In  Figure  3-7  there  is  no  steady  state  solution  because  if  the  steady  state 
values  of  5  and  6  are  placed  into  the  Z  equation,  Z  is  non-zero  below  the  saturation 
Froude  number 

A  parameter  introduced  in  Chapter  II  was  a,  which  allowed  us  to  go  from  a 
multiple  input  system  to  a  single  input  system.  This  parameter  has  a  significant  effect  on 
the  location  of  the  critical  speed  and/or  Froude  number  Figure  3-8  shows  the  relationship 
been  the  critical  speed  ucr  and  the  metacentric  height  zGB  for  different  values  of  oc 

These  three  curves  shown  in  Figure  3-8  can  be  converted  into  a  single  curve  by 
plotting  a  versus  Froude  number,  and  is  shown  in  Figure  3-9  As  can  be  seen  in  this  plot 
the  critical  Froude  number  varies  from  0  81  to  1  22  over  the  allowable  range  of  a 

In  addition  to  the  effect  that  a  has  on  the  location  of  the  critical  Froude  number, 
it  also  effects  the  magnitude  of  the  steady  state  values  at  a  given  Froude  number  and  the 
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STEADY  STATE  VALUES  OF  THETA  vs.  FROUDE  NUMBER 
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Figure  3-5.  Steady  State  Values  of  6  vs  Froude  Number. 
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STEADY  STATE  VALUES  OF  DELTA  vs.  FROUDE  NUMBER 
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Figure  3-6.  Steady  State  Values  of  8  vs  Froude  Number. 
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STEADY  STATE  VALUES  OF  Z  vs.  FROUDE  NUMBER 
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Figure  3-7.  Steady  State  Values  of  Z  vs  Froude  Number. 
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CRITICAL  SPEED  vs.  METACENTRIC  HEIGHT 
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Figure  3-8.  Relationship  Between  Critical  Speed  and  Metacentric  Height  For 

Various  Values  of  a. 
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CRITICAL  FROUDE  NUMBER  vs.  ALPHA 
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Figure  3-9.  Froude  Number  as  a  Function  of  a. 
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location  at  which  the  control  surfaces  saturate  When  comparisons  are  made,  it  can  be 
shown  that  as  the  value  of  a  becomes  more  positive  the  location  of  the  critical  Froude 
number  and  steady  state  values  of  6  and  Z  become  greater  Comparing  the  changes  in  the 
5  steady  state  solutions,  the  shape  of  the  5  pitchfork  does  not  change  however,  the 
location  of  the  critical  point  moves  in  the  direction  of  increasing  Froude  number  as  a 
increases  The  results  of  this  discussion  are  shown  graphically  in  Figures  3-10  through 
3-12. 

An  interesting  result  is  displayed  in  Figure  3-13.  This  figure  depicts  the 
relationship  between  a  and  the  critical  or  saturation  Froude  number.  At  the  lower  Froude 
numbers,  there  is  very  little  difference  between  the  critical  Froude  number  and  the 
saturation  Froude  number  This  demonstrates  that  upon  reaching  a  low  critical  Froude 
number  the  control  surface  immediately  saturates  attempting  to  keep  the  vehicle  stable.  At 
the  higher  Froude  numbers  there  is  a  significant  difference  between  the  critical  Froude 
number  and  the  saturation  Froude  number  This  is  because  as  the  Froude  number 
approaches  2.78,  the  value  used  in  the  control  law  design  conducted  in  Chapter  II,  the 
control  surfaces  will  not  have  to  saturate  to  maintain  stability 
3.      Controllability  Analysis 

In  general,  every  system  of  the  form, 

*:£+*«.  (327) 

can  be  divided  through  a  series  of  transformations  into  four  subsystems: 

1 )  A  controllable  and  observable  part 

2)  An  uncontrollable  and  observable  part. 

3)  A  controllable  and  unobservable  part. 

4)  An  uncontrollable  and  unobservable  part 
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STEADY  STATE  VALUES  OF  THETA  vs.  FROUDE  NUMBER 
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Figure  3-10.  Comparison  of  Steady  State  Values  of  6  for  Different  Values  of  a. 
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STEADY  STATE  VALUES  OF  DELTA  vs.  FROUDE  NUMBER 
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Figure  3-11.  Comparison  of  Steady  State  Values  of  5  for  Different  Values  of  a. 
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FROUDE  NUMBER  vs.  ALPHA 
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This  is  known  as  Kalman's  decomposition  theorem.  Now,  since  the  transfer  function  of 
any  system  is  determined  by  the  controllable  and  observable  subsystems  only,  then  the 
transfer  function  may  contain  less  information  than  is  actually  needed  to  model  the 
complete  system  The  precise  definition  of  controllability  is: 

A  system  is  said  to  be  controllable  if  any  initial  state  x(t0)  can  be  driven  to  any 
final  state  x(tA  using  possibly  unbounded  control  u(t)  in  finite  time  t0<  t  <  tj-. 
From  the  state  equations  3.27,  this  should  depend  only  on  A  and  B,  where  A  is  a  n  x  n 
square  matrix   The  criterion  for  controllability  is  as  follows:  Compute  the  controllability 
matrix 

c  =  [B,AB,A2B A"-XB],  (3.28) 

and  the  system  is  controllable  if  and  only  if  the  rank  of  c  (the  number  of  linearly 
independent  rows  or  columns)  is  n.  Roughly  speaking,  c  shows  how  possible  it  is  to 
change  the  state  of  a  system  using  the  input.  For  a  single  input  system  B  is  n  x  1  and  c  is  a 
square  matrix  The  test  is  then  that  c  be  nonsingular 

det(c)*0.  (3.29) 

Now,  if  the  controllability  matrix  is  formed  using  the  A  and  B  matrices  from 
equation  2.59  then  the  following  matrix  results, 


where. 


c  = 


c„ 

C\2 

c» 

^14 

*21 

c22 

CZ3 

C24 

c31 

C)2 

^33 

^34 

C4, 

CA2 

C43 

C44 

(3.30) 


Cn=0, 


cn  =b2u '  , 


(331) 
(3.32) 
(3.33) 
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c14  =  b{u~ (aua2]u2  +  a2]a22u2 )  +  b2u2 (a23zGB  +  a2]auit2  +a22'u2),  (3  34) 

c2l-b:u2,  (3.35) 

c22  =auuib]+aru}b2,  (336) 

c23  =  b]u'(au'u2  +a2]a]2u2  )  +  b2u2(anzGB  +auauu~  +a]2a22u2 ),  (3  37) 

c,4  =b]u2(auu(au2u2  +a]2a2]u2 )  +  a2]u(anzGB  +  aua]2u2  +a]2a22u2 ))  + 

b2u2  (anzGBauu  +  ana2izGBu  +  a]2u(au2u2  +aua2}u2)+  (3.38) 

a22u(a23zGB  +ana2lu2  +a222u2 )), 

c„=b2u\  (3  39) 

cr  =a^uib]+a„uib2,  (3  40) 

c33  =b]u2(aua^u~  +  a2]a22u2 )  +  b2u2  (a23zaB  +  a2]al2u2  +a22u2 ),  (341) 

c34  =  b]u2(auu(aua2]u2  +ana22u2)+a2lu(a23zGB+al2a2]u2 +a222u2))  + 

b2u~  (anzGBa2Xu  +  a2Za23zGBu  +  ar_u(aua2]u2  +ana22u2)+  (3. 42) 

a22u(a23zGB  +a]2a2]u2  +  a22u2 )), 

c41=0,  (3.43) 

c42=Z>,w:,  (3  44) 

c43  =arni/3A,  +(anu-u)b2u2,  (3  45) 

and  c^  =  by~ (au2u2  +  a2lu(al2u -  u))  +  b2u2 (a]3zGB  +  aua]2u2  +  a22u(a]2u  - u)) .        (3  46) 

Now,  if  the  determinant  of  the  controllability  matrix  is  computed  as  a  function  of 

Froude  number,  then  it  can  be  shown  that  the  system  becomes  uncontrollable  at  the  critical 

Froude  number  The  graphical  results  of  this  computation  are  shown  in  Figure  3-14  for 

different  values  of  the  parameter  a. 

With  the  critical  parameter  that  causes  the  vehicle  to  lose  stability  identified,  i.e., 

the  Froude  number,  and  the  effect  that  the  parameter  a  has  on  the  location  of  this  critical 

value  understood,  the  simplifications  discussed  at  the  beginning  of  this  chapter  can  be 

removed  from  the  E.O.M.  and  an  analysis  on  the  effects  of  incorporating  these  additional 

parameters  can  be  conducted  This  detailed  analysis  will  be  the  topic  of  Chapter  IV. 
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CONTROLLABILITY  vs.  FROUDE  NUMBER 
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IV.  BIFURCATION  ANALYSIS 


A.     ASYMMETRIC  PITCHFORK 

Now  that  the  critical  parameter,  the  Froude  number,  has  been  identified,  the 
simplifications  applied  to  the  E.O.M.  in  the  previous  chapter  can  be  removed  The 
resulting  system  of  equations  is  given  by  equations  2  38  through  2  43  To  determine  the 
existence  of  multiple  steady  state  solutions  for  this  new  system  of  equations  it  is,  once 
again,  necessary  to  perform  a  steady  state  analysis  as  conducted  in  Section  B.2  of  the 
previous  chapter  This  analysis  will  reduce  the  system  of  equations  into  the  following; 

Zwuw--pCDAw\w\  +  (W-B)cosQ  +  Z&u2b  =  0  (4.1) 

Am 

Mj/w--pCDxAAw\w\-(xGW-xBB)cosQ-(zGW-zBB)smQ  +  M&u2b  =  0,  (4.2) 

with  equation  3.13  also  obtained  from  the  steady  state  Z  equation  By  multiplying  equation 
4.1  by  Ms  and  equation  4.2  by  Z6  and  setting  the  resulting  equations  equal  to  each  other 

the  following  equation  is  obtained; 

(ZwMh-MwZh)uw-(-pACD)(M&-xAZh)w\w\  + 

((W  -  B)Mh  +  (xGW  -  xBB)ZJcosQ  +  Z&(zaW-zBB)sinQ  =  0 
This  equation,  along  with, 

w  =  i4tanQ,  (4.4) 

from  equation  3  13  is  the  exact  steady  state  solution  set  with  all  parameters  included 
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Figure  4-1  shows  the  results  obtained  when  equation  4  3  is  solved  numerically,  using  the 
FORTRAN  program  in  Appendix  A,  for  a  small  value  of  xGB  (the  difference  in  the  location 
of  the  center  of  gravity  and  buoyancy),  with  the  submarine  neutrally  buoyant  {W  =  B). 
This  small  value  of  xGB  perturbs  the  pitchfork,  as  seen,  from  the  symmetric  pitchfork 

displayed  earlier  in  Figure  3-5  This  comparison  shows  that  the  critical  or  bifurcation  point 
has  moved  to  a  lower  Froude  number,  with  the  stable  solutions  represented  by  the  solid 
lines  and  the  unstable  solutions  represented  by  the  dashed  lines  in  both  cases.  The  stable 
steady  state  solution,  that  the  vehicle  will  acquire,  is  dependent  on  the  initial  conditions 
that  are  given  to  the  vehicle.  Using  the  results  displayed  in  Figure  4-1  as  an  indication  that 
xGB  may  also  be  a  critical  parameter,  the  subsequent  task  is  to  determine  the  relationship 
between  the  critical  Froude  number  and  xGB  or  f(Fn,xGB). 

B.     BIFURCATION  GRAPHS 

The  first  order  of  business  in  the  attempt  to  identify  f(Fn,xGB)  is  to  simplify  the 

exact  pitchfork  equation  This  can  be  accomplished  by  replacing  the  cosQ  and  sin 6  terms 

in  equation  4  3  with  the  following  Taylor  series  expansions 

w     1  w3      2wu2-w3 

s,nQ  = -— = — (4.5) 

u     2  u  2u 

i         2  o     2  2 

.     .     1  w       2w  -w 
cosQ  =  \--—=  ,      ,  (4.6) 

2  u~  2u~ 

which  will  allow  the  exact  solution  to  be  represented  as  follows, 

Zh(zGW -zBB)w'  +  2u'CDA( Mb-xAZJw\w\ 
-(2uUZ„Mh-MJh)  +  2uzZJzGW-zBB))w-2u-((W-B)Mb,  (4.7) 

+  (xgW-xbB)Z6  +  ((W-B)M&  +  (xgW-xbB)ZJh2=0 
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Figure  4-1.  Exact  Solution  Set  for  Xgb/L=-0.0001. 
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where  CD  now  equals  the  terms  —  pCD  present  in  equation  4.3    Defining  the  following 

parameters 

5,  =  B-W,  (4.8) 

*GB  ~  *G  ~  *B  '  W-") 

ZGB  ~ '  ZG~ '  "fl '  V^    *  " ) 

zGW-zBB  =  zGBW-zBbB  =  5f,  (4. 1 1) 

xGW-xBB  =  xGBW-xBSB=5x,  (4.12) 

X  =  l    Fn\  (4.13) 
and  substituting  them  into  equation  4  7  will  allow  the  coefficients  in  front  of  the  various 
powers  of  w  to  be  written  as  follows; 

u'3:  A4=  Zh(zGBW-zBhB)  =  Z&Cku-W-zBbB),  (4.14) 

w\w\:  A,  =  2u'CDA(Mh-xAZh),  (4.15) 

w:  ^  =-2wV«VZHM5-MRZJ  +  Z6r?iW:^-rs5jA  (416) 

w°:  ^  =-2«V-6sH-Z5rxCT^-xA^  (417) 

w::4)=r-5sM5-Z6rxGfl^-xA»  (4.18) 
The  resulting  form  of  equation  4.7  can  be  written 

A4w^  +  A^w\w\  + A,w  +  A]  +  A^2  =  0,  (4  19) 

which  is  a  generic  pitchfork  equation.  If  it  is  assumed  that  the  submarine  is  neutrally 
buoyant,  and  equation  4.19  is  divided  by  ZhzGBW ,  then  the  above  coefficients  can  be 

redefined  as; 

44  =  1.  (4  20) 

A,=2uCDA(Mh-xAZh)    g      -   ,  (4.21) 

Z*W  &GB 
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j,m_2t(jL*(LMi=M&Ui)t  (4.22) 


(423) 


SZGB 

ZhW 

A  = 

u2 

gZGB 

2g, 

A  = 

u 

2 

XGB% 

gZGB       U' 

Letting  the  critical  parameter  X  be  defined  as 


(4.24) 


*.  =  — ,  (4.25) 


gZGB 

will  allow  equation  4. 19  to  be  put  into  the  following  pitchfork  solution  form; 

w*  +  yhv\w\-2u2(\  +  Q)w  +  2$u2\  +  &kw2  =  0,  (4.26) 

where, 

p  =  ~3aa£  (4.27) 

?=|ft^W_  (428) 

b  zhw 

y  =  2uCDA(Mh-xAZ&)-f-.  (4  29) 

To  determine  if  the  above  manipulation  and  definitions  are  correct  and  valid,  take  the 
symmetric  case  that  was  developed  in  the  previous  chapter,  i.e.,  xGB  =  0.  This  causes 

equation  4.26  to  reduce  to 

w(w2  +  yX\w\  -  2u2(\  +  C,X)J  =  0,  (4.30) 

resulting  in  a  solution  always  occurring  at  w=0,  with  two  more  solutions  appearing  at  the 
bifurcation  value 

1  +  C\,  =  0.  (4.31) 

Rearranging  and  substituting  the  appropriate  values  into  equation  4.31  will  yield 
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(4.32) 


which  is  identical  to  equations  3  9  and  3.20,  thus  verifying  the  manipulation. 

To  find  the  relationship  between  Fn  and  xGB  the  function  h(u,\,$),  and  its  partial 
derivative  hH(u\\,$),  must  be  determined  These  two  functions  must  then  be  set  equal  to 
zero  and  to  each  other  in  order  to  obtain  a  function  independent  of  w.  As  seen,  equation 
4.30  is  much  too  complex  to  complete  these  requirements,  but  by  neglecting  the  terms 
Xm  |h  [  and  \w~  in  equation  4  26  a  simplified  pitchfork  of  the  form, 

wi-2uz(\  +  Q)w  +  2$u2\  =  0,  (4.33) 

or  in  general  terms  h(w,\,$)  =  0,  can  be  obtained  which  can  be  manipulated  to  meet  the 
above  mentioned  requirements,  thereby  determining  the  critical  (\,$)  curve.  By  taking 
the  partial  derivative,  /7M„  of  equation  4.33,  the  following  equation  results, 

hw  =  3w2-2ul(l  +  &).  (4.34) 

This  equation  can  be  set  equal  to  zero,  yielding 

wa=!«Yl  +  SU  (4.35) 

which  can  be  substituted  into  equation  4.33  to  obtain  the  value  of  w  where  the  function 
and  it's  partial  derivative  is  equal  to  zero; 

w  = (4  36) 

This  value  of  w  can  now  be  substituted  back  into  equation  4.35  to  obtain  the  desired 
/(VPJ.or 

27$2X2=&u2(\  +  {,\f.  (4-37) 

Attempting  to  solve  this  equation  analytically  for  X  as  a  function  of  (3,  to  determine 
the  shape  of  this  curve  is  difficult  Therefore,  an  approximation  to  this  equation  is  needed 
for  small  values  of  (3   By  rearranging  equation  4.37  into  the  form 
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y>W,  (4  38) 

an  equation  that  determines  (3  as  a  function  of  X  is  now  produced.  If  this  function  is 
represented  in  a  Taylor  series  expansion  and  simplified,  then  equation  4.38  can  be  written 
as 

x=x0+|rp2/r«2c4;/3,  (4.39) 

which  can  be  considered  as  an  approximation  to  the  (X,$)  curve,  with  equation  4.37  as 
the  exact  (X,$)  curve  Through  inspection  it  can  be  seen  that  the  general  shape  of  this 
curve  is  that  of  a  cusp. 

Figure  4-2  shows  a  comparison  of  four  bifurcation  curves.  The  curves  represent  the 
relationship  between  the  critical  speed  U,  which  can  be  converted  into  the  non  dimensional 
parameter  X  using  equation  4.25,  and  the  value  of  xGB,  which  also  can  be  converted  into  a 
non  dimensional  parameter  (3  using  equation  4.27  The  curves  labeled  exact  bifurcation 
set,  pitchfork  bifurcation  set,  exact  cusp  and  approximate  cusp  come  from  the  numerical 
solution  to  equations  4.7,  4.26,  4.37  and  4.39  respectively.  Observe  that  in  the  general 
vicinity  of  xGB=Q  and  near  the  critical  speed  the  shape  of  each  of  these  curves  is  indeed 
that  of  a  cusp,  with  the  peak  occurring  at  the  critical  speed  of  1  9  fps  Now,  with  the 
comparison  conducted  in  Figure  4-2  complete,  equation  4.7,  the  exact  pitchfork  solution 
will  be  used  for  the  remainder  of  this  analysis. 

Returning  to  equation  4.7,  it  can  be  seen  that  there  is  only  one  additional  parameter, 
CD,  that  the  cusp  curve  produced  using  equation  4  39  does  not  incorporate  (because  zGB 
and  xGB  are  incorporated  into  the  non  dimensional  parameters  X  and  (3)  It  is  therefore 
necessary  to  conduct  a  sensitivity  analysis  to  determine  the  effect  that  CD  has  on  the  shape 
of  the  cusp  curve.  Figure  4-3  displays  the  results  obtained  when  this  sensitivity  analysis 
was  conducted 
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Bifurcation  Sets  (CD=0,  Zgb=0.1) 
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Figure  4-2.  Comparison  of  Bifurcation  Curves. 
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Bifurcation  Sets  (Zgb=0.1) 
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Figure  4-3.  Comparison  of  Exact  Bifurcation  Set  for  Different  Values  of  Drag 

Coefficient. 
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As  seen  in  the  plot,  the  drag  coefficient  has  no  effect  on  the  location  of  the  peak  of 
the  cusp  The  drag  coefficient,  however,  does  modify  the  slope  of  the  cusp  as  displayed  by 
the  four  cases  shown  in  Figure  4-3. 

C.     SOLUTION  SETS 

In  this  section,  the  equilibrium  solutions  to  equations  4.3  and  4  5  will  be  investigated. 
As  determined  in  the  previous  section,  the  critical  parameters  that  control  the  location  and 
shape  of  the  pitchfork  solution  sets  are  X,  P  and  the  drag  coefficient  CD.  By  running  the 
FORTRAN  program  of  Appendix  A  for  several  values  of  xGB  L  the  shape  and  trend  of 
the  solution  sets  can  be  determined  The  result  of  this  comparison,  for  a  CD  =  0.0,  is 
displayed  in  Figure  4-4  This  plot  supports  the  symmetry  about  xGB  =0  displayed  in  the 
earlier  cusp  curves  Observe  that  the  critical  or  bifurcation  point  moves  to  a  lower  value  of 
X  as  xGB  increases,  and  that  for  low  values  of  X  the  equilibrium  solutions  converge  to  the 
same  value  independent  of  the  value  ofxGB. 

When  the  comparison  is  made  regarding  the  effect  that  the  drag  coefficient  has  on  the 
solution  sets,  for  a  given  value  of  xGB,  Figure  4-5  is  produced  This  variation  in  drag 
coefficient  for  a  constant  value  of  xGB  has  a  similar  effect  on  the  movement  of  the 
bifurcation  point  as  fixing  the  drag  coefficient  and  varying  xGB.  The  difference  is,  however, 
evident  by  studying  the  stable  solutions  The  stable  solutions  for  non-zero  drag 
coefficients  become  more  linear  as  the  value  of  CD  increases  This  results  from  the 
increased  effect  of  the  w\w\  term  in  equation  4  3  as  the  value  of  CD  increases 

Although  Figure  4-5  shows  the  actual  solution  to  the  exact  pitchfork  equation  it  is 
not  realistic  since  the  effect  of  dive  plane  saturation  is  not  considered  Once  the  dive 
planes  saturate,  the  equations  used  to  obtain  the  steady  state  solutions  displayed  in 
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Solution  Set  (CD  =  0.0,  Zgb=0.1) 
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Figure  4-4.  Comparison  of  Exact  Solution  Set  for  Different  Values  of  Xgb/L. 
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Solution  Set  (Zgb=0.1,  Xgb/L= -0.0001) 
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Figure  4-5.  Comparison  of  Exact  Solution  Set  for  Different  Values  of  CD. 
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Figures  4-4  and  4-5  are  no  longer  valid  Instead,  equations  4  1  and  4  2  are  modified,  still 
remembering  that  the  assumption  of  neutral  buoyancy  remains  in  affect,  resulting  in  the 

following  two  equations, 

Zjuw  -  CDAw\w\  +  u-Zbhsat  =  0,  (4.40) 

Mjtw  -  CDxAAw\*>\  -  WxGB  cosQ-  zGBWsinQ  +  M&u'bsal  =  0 .  (4.41 ) 

To  obtain  an  equation  that  can  be  used  to  compute  the  steady  state  0  solution  while 
taking  into  consideration  dive  plane  saturation  requires  solving  equation  4.40  for  w  and 
then  substituting  this  value  into  equation  4  41  which  should  yield  an  equation  independent 
of  w.  The  ability  to  perform  this  manipulation  is  hindered  by  the  inclusion  of  the  absolute 
value  term  in  both  of  the  above  equations.  However,  if  the  absolute  value  terms  are 
neglected  by  considering  the  case  of  CD=0.0,  then  a  single  equation  for  the  steady  state 
value  of  6  may  be  obtained 

Neglecting  the  drag  affects,  and  performing  the  algebra  will  yield  an  equation  of  the 
following  form, 

(  MwZb  -  MbZJu%al  +  ZW(xGB  cosQ  +  zGB  sin  Q)  =  0,  (4.42) 

which  may  be  used  to  determine  the  values  of  0  once  the  dive  planes  have  saturated.  The 
validity  of  neglecting  the  drag  terms  is  questionable,  and  an  analysis  that  incorporates  CD 
will  be  conducted  later  in  this  section  once  the  effect  of  varying  xGB  is  determined  Figure 
4-6  displays  the  changes  that  occur  to  the  solution  set  of  Figure  4-1  when  saturation  is 
considered.  Notice  that  saturation  occurs  when  the  solution  set  is  at  approximately  nine 
degrees,  and  that  when  multiple  solutions  occur,  the  stable  equilibrium  states  are 
represented  by  the  saturation  portions  of  the  solution  set 
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Solution  Set  (CD=0.0,  Zgb=0.1,  Xgb/L= -0.0001) 
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Figure  4-6.  Exact  Solution  Set  for  Xgb/L=-0.0001  With  Saturation  Included. 
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When  the  saturation  condition  is  applied  to  the  previously  displayed  solution  sets  for 
different  values  of  xGB  (Figure  4-4),  the  results  shown  in  Figure  4-7  are  produced  Items 
of  interest  in  Figure  4-7  are  that: 

1)  The  bifurcation  point  on  each  solution  set  moves  to  a  lower  value  of  X  as 
the  absolute  value  of  xGB  increases 

2)  The  range  where  multiple  solutions  exist  decreases  as  the  absolute  value  of 
xGB  increases. 

3)  The  amount  of  trim  produced,  in  still  water,  by  the  given  value  of  xGB  can  be 
obtained  from  the  point  where  the  upper  an  lower  solution  branches 
converge 

Incorporating  the  saturation  condition  into  the  solution  sets  raises  some  additional 
concerns  To  begin,  what  if  the  point  of  dive  plane  saturation  occurs  after  the  bifurcation 
points  identified  in  Figure  4-49  If  this  condition  did  occur,  referring  to  Figure  4-1  which 
displays  the  stable  and  unstable  solution  branches,  there  would  be  at  least  one  value  of  X 
where  there  were  three  stable  solution  branches.  To  determine  if  this  condition  is  possible 
a  comparison  of  solution  sets  with  and  without  saturation  is  necessary  This  comparison  is 
shown  in  Figure  4-8  As  this  plot  depicts,  the  range  between  the  bifurcation  point  and 
saturation  point  increases  as  the  value  of  xGB  increases.  This  result  will  allow  the  fore 
mentioned  concern,  of  three  stable  solution  branches,  to  be  disregarded. 

A  second  concern  is  the  validity  of  the  previously  developed  bifurcation  cusp  curves 
Since  the  exact  bifurcation  equation,  equation  4-7,  does  not  include  saturation,  the  cusp 
curves  produced  using  this  equation  will  not  accurately  predict  the  occurrence  of  multiple 
solutions  However,  if  equation  4.42  is  expanded  in  a  Taylor  series,  it  can  be  manipulated 
into  a  form  that  will  allow  a  cusp  curve,  which  includes  saturation,  to  be  developed. 
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Figure  4-7.  Comparison  of  Exact  Solution  Set  for  Different  Values  of  Xgb/L  With 

Saturation  Included. 


66 


20 


15- 


10 


5- 


u 
u 

— 

u 

s 


-51- 


■10  h 


■15  ^ 


.2 


zl 


Solution  Set  (CD =0.0,  Zgb=0.1) 


4  0.6  0.8 


for  Dive  Plane  Saturation 

-  Xgb/L  =0.0001 

-  Xgb/L  =0.0003 
-.  Xgb/L  =0.0005 


1  1.2  1.4  1.6  1.8  2 

SQRT(LAMBDA) 

Figure  4-8.  Comparison  of  Bifurcation  Points  and  Dive  Plane  Saturation  Points. 
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Performing  this  manipulation  will  produce  the  following  equation; 

-x0BZwW  +  zaBW2^,  M4J 

(ZiMw-MtZJ(SM(l+(Zs/Zj5^D 

By  converting  u:  and  xGB  into  the  non  dimensional  parameters  X  and  (i  the  cusp 

curve  shown  in  Figure  4-9  can  be  produced.  This  cusp  displays  the  saturation  point  for 

both  solution  branches,  and  as  seen  it  differs  quite  extensively  from  the  exact  bifurcation 

cusp  without  saturation  If  Figure  4-9  is  replotted  to  allow  the  peak  area  of  the  cusp  to  be 

more  accurately  represented,  then  Figure  4-10  is  produced.  By  using  these  two  cusp 

curves  then  the  general  form  of  the  resulting  solution  set  for  any  path  through  the  cusp  can 

be  predicted  It  is  this  prediction  that  the  subsequent  section  will  address 

D.     PATH  FORMULATION 

The  ability  to  accurately  predict  how  the  vessel  will  respond  to  various  changes  in 
operating  conditions  is  critical  for  proper  control  of  the  vessel  It  is  this  prediction  that 
necessitates  the  need  for  path  analysis  to  be  conducted  If  the  ballast  condition  on  the 
vessel  is  held  constant,  xGB  =-0.0001,  and  the  speed  of  the  vessel  is  varied,  path  A  of 
Figure  4-11  is  obtained.  Neglecting  for  the  moment  the  saturation  cusp  shown  in 
Figure  4-11,  it  would  be  expected  that  a  single  steady  state  solution  would  exist  for  all 
values,  of  the  parameter  sqrt(X),  greater  than  the  point  where  the  path  intersects  the  solid 
cusp  curve  At  the  speed  where  the  path  intersects  the  cusp  the  bifurcation  point  occurs, 
and  for  speed  values  which  place  the  path  inside  the  cusp  multiple  steady  state  solutions 
occur  The  results  of  this  path  analysis  are  depicted  in  Figure  4-1,  presented  earlier  in  this 
chapter. 

If  attention  is  returned  to  the  path  through  the  saturation  cusp  of  Figure  4-11,  then 
the  same  type  of  predictions  made  for  the  path  through  the  solid  cusp  can  be  made  for  the 
this  cusp,  with  some  slight  differences  As  the  speed  decreases  the  path  intersects  the 
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Figure  4-9.  Bifurcation  Cusp  With  Dive  Plane  Saturation  Included. 
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Bifurcation  Sets  (CD  =  0.0) 
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Bifurcation  Sets  Path  Formulation 


Figure  4-11.  Path  Formulation  for  Xgb/L=-0.0001. 
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upper  portion  of  the  dotted  cusp  at  an  approximate  value  of  1.08  At  this  value  saturation 
occurs  on  one  of  the  solution  branches  but  there  remains  only  one  steady  state  solution 
since  the  path  has  not  yet  entered  the  cusp  As  the  speed  continues  to  decrease,  the  path 
intersects  and  enters  the  saturation  cusp  at  an  approximate  value  of  0.99  with  saturation  of 
the  other  solution  branch  and  the  existence  of  multiple  steady  state  solutions  occurring. 
These  results  are  evident  in  Figure  4-6 

Now,  if  the  speed  of  the  vessel  is  held  constant,  and  the  loading  conditions  on  the  hull 
are  changed,  then  the  paths  shown  in  Figure  4-12  are  obtained  Path  £,  at  a  Froude 
Number  of  1.1,  is  outside  the  cusp,  therefore  only  a  single  steady  state  solution  for  each 
value  of  xGB  L  should  exist  While  path  C,  at  a  Froude  Number  of  10,  passes  through  the 
cusp  at  ±  0  0001  values  oixGBL  resulting  in  single  steady  state  solutions  for  values  less 
than  or  greater  than  ±  0.0001,  while  values  between  ±  0  0001  result  in  multiple  steady 
state  solutions  Through  numerical  computation  the  solution  set  as  a  function  of  xGBL  is 
obtained  with  the  results  shown  in  Figure  4-13. 

As  seen  in  Figure  4-13,  path  B,  represented  by  the  dotted  line,  does  in  fact  result  in  a 
single  steady  state  solution  throughout  the  range  ofxGB/L,  however,  for  path  C  this  is  not 
the  case.  Path  C,  represented  by  the  solid  line,  exhibits  multiple  solutions  between  the 
xGB  L  values  of  ±  0  0001,  as  predicted  This  is  a  hysteresis  curve,  with  the  unstable 
solutions  occurring  between  the  corners  of  the  curve  If  the  value  of  xGB/L  began 
at  -0.0005  and  progressed  to  the  positive  value  of  0.0001,  the  solution  set  would  remain 
on  the  curve  with  6  taking  on  the  values  shown  Once  the  value  ofxGBL  becomes  greater 
than  0.0001  the  solution  jumps  from  -9°  to  7.5°.  This  same  jump  would  be  evident  if  the 
value  of  xGB  L  began  positive  and  progressed  negative  These  jumps  indicate  that  if  the 
dive  planes  are  held  constant,  then  the  vessel  will  instantaneously  change  its  pitch  angle 
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Bifurcation  Sets  Path  Formulation 
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Figure  4-12.  Path  Formulation  at  a  Constant  Speed  With  Xgb/L  Varying. 
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Solution  Set  for  Path  Formulation 
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To  counter  this  effect  the  dive  plane  angle  must  be  reversed  It  is  this  dive  plane  reversal 
that  the  following  chapter  will  discuss. 

Before  proceeding  to  the  next  chapter  it  is  necessary  to  determine  the  effect  that  drag 
has  on  the  previously  developed  bifurcation  cusps  Using  equations  4 .40  and  4.41,  and 
substituting  equation  4.4  into  both  equations  will  result  in  the  following  two  equations, 

Zhm2  tan  0-C^w2  tan  6|tan  6|  +  Z^S,,,  =0,  (4.44) 

M„uw-CDxAAu2  tan  6|tan6J-  WxGB  cos  Q-  zGBW  sin  0+  Msu28sal  =  0      (4.45) 

By  using  the  MATLAB  function  /zeros  contained  in  the  program  of  Appendix  A,  the 
value  of  6  that  solves  equation  4.44  can  be  determined  as  a  function  of  speed  This 
relationship  can  then  be  used  to  solve  equation  4  45  to  obtain  the  relationship  between  xGB 
and  u 

With  this  relationship  obtained  the  cusp  curves  of  Figure  4-14  can  be  plotted  It  can 
be  seen  that  for  a  given  value  of  xGB/L  increasing  the  drag  coefficient  forces  the 
bifurcation  point  to  a  lower  speed  An  interesting  item  to  note  is  the  diamond  shape  peak 
on  these  cusps.  This  is  a  result  of  saturation  not  occurring,  for  a  small  range  of  speeds,  in 
the  small  region  around  xGB/L=0.  Path  formulation  through  these  cusps  can  be  developed 
as  shown  earlier  in  this  section. 
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Figure  4-14.  Bifurcation  Cusp  With  Saturation  Included  for  Different  Values  of  CD. 
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V.  DIVE  PLANE  REVERSAL 


A.  STERN  PLANE  REVERSAL 

The  ability  of  a  submarine  to  avoid  detection  and  perform  its  mission  is  dependent  on 
its  ability  to  maintain  ordered  depth.  Modern  asymmetric  submarines  exhibit  an  inherent 
phenomenon  known  as  dive  plane  reversal,  which  is  difficult  to  deal  with.  To  introduce 
this  physical  occurrence,  refer  to  Figure  5-1  [Ref.  14]  and  consider  the  case  of  wanting  to 
dive  the  vehicle  to  a  deeper  depth.  At  moderate  to  high  speeds  the  stern  planes  would  be 
deflected  by  an  angle  5S,  this  would  in  turn  cause  a  force  and  moment  to  be  developed  due 
to  the  angle  of  attack  of  the  planes.  The  vessel  will  respond  by  pitching  downward,  and 
will  assume  some  angle  of  attack  to  the  fluid  flow  causing  hull  forces  and  moments  to 
develop  The  stern  plane  and  hull  forces  bring  the  ship  to  some  steady  state  pitch  angle 
which  is  obtained  when  the  restoring  moment  and  pitching  moments  balance  The  vessel 
now  has  the  ability  \ofly  itself  to  the  ordered  depth  with  the  normal  hull  force  assisting  in 
pushing  the  vessel  to  a  deeper  depth. 

At  low  speed  we  would  expect  the  same  to  be  true,  however  it  is  not.  If  the  stern 
planes  are  set  to  an  angle  5S ,  this  would  cause  a  force  and  moment  to  be  developed  due  to 
the  angle  of  attack  of  the  planes,  the  vessel  will  pitch  downward  creating  a  pressure  field 
around  the  vessel  due  to  the  angle  of  attack  with  the  fluid  resulting  in  hull  forces  and 
moments  to  develop.  The  vessel  once  again  steadies  out  at  some  steady  state  pitch  angle 
when  the  restoring  moment  and  pitching  moments  balance.  Since  the  speed  of  the  vessel  is 
not  as  great  as  in  the  previous  case,  the  hull  forces  and  moments  are  not  as  large,  therefore 
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Figure  5-1.  Stern  plane  Reversal 
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instead  of  the  vessel  flying  itself  to  the  ordered  depth  with  the  hull  force  assisting,  the 
stern  plane  force  actual  pushes  the  vessel  upward  causing  it  to  rise  instead  of  dive,  with 
some  non-zero  pitch  angle.  [Ref  14]  This  condition  was  evident  in  the  simulation  shown 
in  Figure  3.3. 

With  this  description  of  the  vehicle  response  complete,  it  is  obvious  that  to  dive  the 
vessel  at  low  speeds,  the  opposite  stern  plane  angle  must  be  ordered.  This  will  result  in  the 
vessel  being  pushed  to  the  ordered  depth  at  some  bow  up  attitude  This  required  shift  in 
dive  plane  angle  is  what  is  known  as  STERN  PLANE  REVERSAL. 

To  tie  together  the  previously  completed  bifurcation  analysis  with  the  physical 
occurrence  just  discussed,  Figure  5-2  is  produced.  Figure  5-2  presents  the  pitch  angle  and 
stern  plane  angle  required  for  straight  and  level  flight  with  the  bow  planes  centered  As  the 
vessel's  speed  decreases,  the  metacentric  moment  (Med)  makes  it  increasingly  difficult  to 

maintain  a  positive  pitch  angle  on  the  hull  The  stern  plane  angle  is  negative  trying  to 
produce  a  pitch-up  moment,  but  the  downward  force  of  the  stern  planes  requires  an  even 
larger  pitch  angle  [Ref.  14]  Between  the  speeds  of  1.08  and  1.18  knots,  the  required  stern 
plane  angle  is  beyond  the  normal  operating  range  and  the  vessel  cannot  be  flown  straight 
and  level  Below  the  critical  speed,  the  stern  plane  and  vessel  pitch  angle  must  reverse  sign 
to  maintain  neutral  trim. 

In  a  final  effort  to  show  the  physical  significance  of  stern  plane  reversal,  the  force 
input  to  the  stern  planes  is  plotted.  If  the  control  gains  are  allowed  to  change  as  the  speed 
of  the  vessel  changes.  Figure  5-3  is  produced  This  plot  shows  the  non-dimensional  force 
per  unit  depth  error  input  to  the  stern  plane  as  a  function  of  speed  It  is  basically  the 
control  system  gain  used  to  set  the  stern  plane  angle  As  can  be  seen,  as  the  critical  speed 
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Figure  5-2.  Variation  of  Pitch  and  Stern  plane  Angles  as  a  Function  of  Speed. 
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Figure  5-3.  Control  System  Force  Input  to  the  Stern  planes. 
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is  approached  the  value  of  input  tends  to  infinity  requiring  the  stern  planes  to  operate 
outside  normal  values  As  the  speed  passes  the  critical  value  the  sign  of  this  input  changes 
requiring  the  stern  planes  to  operate  in  the  opposite  direction.  This  reversal  in  sign  is 
evident  of  the  stern  plane  reversal  phenomenon  shown  and  discussed  earlier  in  this 
chapter. 

B.     BOW  PLANE  REVTRSAL 

The  reversal  effect  that  was  shown  to  occur  with  the  stern  planes  can  also  be  shown 
to  occur  with  the  bow  planes.  If  the  value  of  a  is  allowed  to  tend  to  oo,  indicating  no  stern 
plane  control  and  only  bow  plane  control,  then  the  same  analysis  conducted  throughout 
this  thesis  for  stern  plane  reversal  can  be  performed  for  bow  plane  reversal.  If  the  critical 
Froude  number,  (sqrt(X)),  is  plotted  as  a  function  of  a,  then  the  results  shown  in 
Figure  5-4  are  produced  This  plot  shows  that  for  bow  plane  control  only,  the  critical 
value  approaches  2.05  as  a  approaches  oo 

To  show  that  the  same  type  of  results  can  be  obtained  for  the  bow  plane  analysis,  the 
case  of  neutral  buoyancy,  xGB  =  0  and  dive  plane  saturation,  with  bow  plane  control  only 
was  investigated  Referring  to  Figure  5-5,  it  can  be  seen  that  at  a  Froude  number  of  2.05 
multiple  steady  state  solutions  occur  These  results  are  similar  to  the  results  displayed  in 
Figure  3-5  The  major  differences  are  that  the  critical  point  occurs  at  a  Froude  number  of 
2.05  vice  1.05,  and  the  maximum  value  of  0  is  only  ±  4.5°  vice  ±  9.5°  The  reason  for 
these  changes  is  due  to  the  values  of  Z5  and  MSt  being  only  one-half  the  respective  value 

for  the  stern  planes 
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Figure  5-4.  Identification  of  Critical  Froude  Number  for  Bow  plane  Control. 
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Figure  5-5.  Steady  State  Values  of  6  for  Bow  plane  Control. 
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C.     BIAS  EFFECTS 

As  a  final  area  of  analysis,  the  condition  of  operating  near  the  surface,  such  as  during 
periscope  operations,  must  be  investigated  to  determine  the  effects  this  will  have  on  the 
previously  developed  bifurcation  cusps  Operating  near  the  surface  will  require  the 
incorporation  of  a  suction  force  and  moment  into  the  equations  of  motion.  This  force  and 
moment  is  a  function  of  the  distance  away  from  the  surface,  and  the  angle  the  vessel 
makes  with  the  surface.  Therefore,  as  a  simple  model  of  suction  effects  consider  the  case 
of  operating  with  excess  buoyancy  To  conduct  the  analysis  requires  the  use  of  equations 
4.1,  4.2  and  4  4  Performing  the  same  substitution  and  solution  technique  as  was  used  in 
Chapter  V  to  obtain  the  saturation  curves  that  incorporated  the  drag  coefficient  will  result 
in  solving  the  following  equations  to  obtain  the  bifurcation  cusps  that  would  quasi-model 
surface  effects; 

Zjiw  --pCpAu2  tan  ^n6\  +  (W  -  B)cos6+  Zsu28sal  =  0,  (5.1) 

Mwuw — pCDxAAu 2  tan  fljtan  0)  -{xGW  -  xBB)cos& 
-{zGW-zBB)sme+Mju285a,  =  0 
If  the  value  of  CD  is  set  equal  to  zero  then  the  saturation  cusp  produced  in 
Figure  4-10  can  be  compared  to  the  results  obtained  for  several  cases  of  excess  buoyancy. 
Since  it  is  the  saturation  cusp  that  controls  the  location  of  the  bifurcation  point,  as  shown 
in  previous  the  chapter,  the  non  saturation  cusp  will  be  neglected  for  the  remainder  of  the 
analysis  Figure  5-6  shows  the  results  obtained  as  excess  buoyancy  is  increased  As  can  be 
seen,  the  cusp  peak  moves  to  down  and  to  the  right  occurring  at  a  lower  speed  and  at  a 
non  zero  value  ofxGB,L.  The  body  of  the  cusp  also  shifts  to  the  right,  as  seen,  with  the 
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Figure  5-6.  Bias  Effects  Caused  by  Operating  Near  the  Surface. 
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entire  cusp  occurring  in  the  positive  region  of  the  xGB/'L  range  for  a  value  of  0.05  % 
excess  buoyancy. 

To  determine  how  the  value  of  drag  coefficient  effects  the  cusp  curves,  consider  the 
case  of  CD  =  0.3  shown  in  Figure  5-7  These  cusps  were  produced  for  the  same  values  of 
excess  buoyancy  used  in  Figure  5-6  Notice  that  the  cusp  curves  shift  to  the  right  by  the 
same  value  seen  in  Figure  5-6,  but  the  peak  value  now  occurs  at  an  even  a  lower  speed 

In  order  to  observe  the  sensitivity  of  the  biased  cusp  to  drag  only,  consider  the  case 
of  setting  the  excess  buoyancy  to  0.01  %  and  varying  the  drag  coefficient  shown 
Figure  5-8  These  results  indicate  that  increased  drag  causes  the  critical  speed  to  occur  at 
a  lower  value  They  also  show  that  drag  tends  to  counter  the  biasing  effect  that  excess 
buoyancy  was  previously  show  to  have  This  is  evident  by  the  peak  point  drifting  to  the 
left  as  the  value  of  drag  coefficient  increases. 

In  summary,  it  can  be  stated  that  operating  near  the  surface  will  tend  to  reduce  the 
critical  speed  determined  in  previous  chapters  and  the  vessel  will  have  to  be  placed  in  a 
different  trim  condition  to  counter  surface  effects,  and  that  increasing  drag  also  decreases 
the  critical  speed  but  tends  to  counter  the  trim  condition  required 
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Figure  5-7.  Bias  Effects  Caused  by  Operating  Near  the  Surface  Considering  Drag. 
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BIASED  BIFURCATION  CUSP  (0.01  %  excess  Bouyancy) 
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Figure  5-8.  Bias  Effects  Caused  by  Drag  While  Operating  Near  the  Surface. 
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VI.  CONCLUSIONS  AND  RECOMMENDATIONS 


A.  CONCLUSIONS 

•  A  method  for  analyzing  submarine  depth  keeping  and  dive  plane  reversal  has  been 
developed  and  presented  This  analysis  identified  the  three  phenomena,  a  positive  real 
eigenvalue,  a  steady  state  pitchfork  bifurcation,  and  an  uncontrollable  system,  that  lead  to 
vessel  instability 

•  The  critical  parameters  have  been  identified  that  control  the  bifurcation.  These 
parameters  combine  speed,  loading  conditions  and  hydrodynamic  coefficients  into 
non-dimensional  values  that  may  be  applied  to  other  models  or  vessels 

•  Since  the  non-dimensional  Froude  numbers  are  based  on  the  metacentric  height 
and  LCG-LCB  separation,  Froude  scaling  can  be  used  to  apply  this  analysis  to  full  scale 
vessels,  allowing  for  the  reduction  of  model  testing  during  design  phases 

B.  RECOMMENDATIONS 

Some  recommendations  for  further  research  are  as  follows: 

•  Develop  and  incorporate  into  the  bifurcation  analysis  a  more  complete  expression 
for  the  free-surface  effects. 

•  Conduct  a  Hopf  bifurcation  analysis  to  identify  the  critical  parameters  that  lead  to 
vehicle  instability  at  high  speeds 
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APPENDIX  A 


COMPUTER  PROGRAMS 

%  THIS  PROGRAM  IS  THESIS  l.M.  IT  SIMULATES  THE  RESPONSE  OF  THE 

DTRC  SUBOFF 

%  MODEL  IT  PLACES  THE  POLES,  CALCULATES  THE  CRITICAL  SPEED,  AND 

USES 

%  PROGRAM  THESISDE  M  TO  SIMULATE  THE  VERTICAL  PLANE  REPONSE 

clg; 

global  zg  m  zgb  u  U  Zq  Zw  Zdlt  Mq  Mw  Bl  Mdlt  Kl  MbmxL  Zval  Mval  rho 

%  ESTABLISH  GEOMETRIC  PARAMETERS 

W=1556.2363,B1  =  1556.2363; 

Iy=561.32; 

g=32.2; 

m=W/g; 

rho=1.94; 

L=13.9792, 

xg=0;zg=0.1, 

zgb=0.1; 

%  BEAM  AND  LENGTH  DEFINITION  FOR  TRAPAZIODAL  RULE  USE 

bm=[0  485  658  778  871  945  1.01  1.06  1  18  1.41  1.57  1.66  1.67  1.67  1.67  1  63 
1.37  919  448  195  .188  168  132  053  0], 

xl=[0  .1  .2  .3  .4  .5  .6  .7  1  2  3  4  7.7143  10  15.1429  16  17  18  19  20  20.1... 
20.2  20.3  20  4  20.4167], 
xl=(L/20)*xl; 
xL=xl-L/2; 

%  FACTORS  TO  CONVERT  FROM  NON-DIMENSIONAL  VALUES  TO 
DIMENSIONAL  VALUES 

ndl=5*rho*LA2, 
nd2=  5*rho*LA3, 
nd3=.5*rho*LA4; 
nd4=5*rho*LA5; 
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AJpha=0; 

%  NON-DIMENSIONAL  HYDRODYNAMIC  COEFFICIENTS 
Zqdnd=-6.33e-4;Zwdnd=-1.4529e-2;Zqnd=7.545e-3;Zwnd=-1.391e-2; 

Zds=(-5.603e-3);Zdb=0.5*(-5.603e-3);Zdltnd=(Zds+Alpha*Zdb); 
Mqdnd=-8  8e-4;Mwdnd=-5 .6 1  e-4;Mqnd=-3  702e-3 ,Mwnd=  1  0324e-2, 
Mds=(-0.002409);Mdb=0.5*(0.002409);Mdltnd=(Mds+Alpha*Mdb); 

%  CONVERTION  OF  COEFFICENTS  INTO  DIMENSIONAL  VALUES 
Zqd=nd3  *Zqdnd,Zwd=nd2*Zwdnd;Zq=nd2*Zqnd;Zw=nd  1  *Zwnd; 
Zdlt=ndl*Zdltnd, 

Mqd=nd4*Mqdnd;Mwd=nd3*Mwdnd;Mq=nd3*Mqnd;Mw=nd2*Mwnd; 
Mdlt=nd2*Mdltnd; 

%  FORMULATION  OF  A  AND  B  MATRIX  ELEMENTS 

Dv=(m-Zwd)*(Iy-Mqd)-(m*xg+Zqd)*(m*xg+Mwd); 

al  lDv=(Iy-Mqd)*Zw+(m*xg+Zqd)*Mw; 

a  1 2Dv=(Iy-Mqd)*(m+Zq)+(m*xg+Zqd)*(Mq-m*xg); 

al3Dv=-(m*xg+Zqd)*W; 

b  1  Dv=(Iy-Mqd)*Zdlt+(m*xg+Zqd)*Mdlt, 

a2 1  Dv=(m-Zwd)*Mw+(m*xg+Mwd)*Zw; 

a22Dv=(m-Zwd)*(Mq-m*xg)+(m*xg+Mwd)*(m+Zq), 

a23Dv=-(m-Zwd)*W; 

b2Dv-(m-Zwd)*Mdlt+(m*xg+Mwd)*Zdlt, 

al  l=al  lDv/Dv,al2=al2Dv/Dv;al3=al3Dv/Dv, 
a2 1  =a2 1  Dv/Dv;a22=a22Dv/Dv,a23=a23Dv/Dv, 
b  1  =b  1  Dv/Dv,b2=b2Dv/Dv; 

%  ESTABLISHMENT  OF  POLE  LOCATIONS 

pl=[-3  -31  -32-33], 

%  CALCULATION  OF  CRITICAL  SPEED 

Ucr=sqrt(-(al3*b2-a23*bl)*zgb/(all*b2-a21*bl)) 

%  INPUT  OF  SPEED  WHICH  THE  VESSEL  WILL  USE  IN  THE  SIMULATION 

U=input('enter  speed  of  vehicle  ') 

%  CALCULATION  OF  FROUDE  NUMBER 

Fn=sqrt(UA2/(zgb*g)) 
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%  CALCULATION  OF  MATRICES  FOR  POLE  PLACEMENT  AND 
CONTROLABILLITY  ANALYSIS 

u=5; 

a=[0  0  1  0;al3*zgb  al  l*u  al2*u  0;a23*zgb  a21*u  a22*u  0; 

-u  10  0]; 

al=[0  0  1  0;al3*zgball*Ucral2*Ucr0;a23*zgba21*Ucra22*Ucr0;... 

-Ucr  1  0  0]; 

b=[0;bl*uA2;b2*uA2;0], 
bl=[0;bl*UcrA2;b2*UcrA2;0], 

Kl=place(a,b,pl), 
C=rank(ctrb(al,bl)); 

H=(Mw*Zdlt-Zw*Mdlt)/(zgb*B  1  *ZdIt); 

thetal=acos(H*UA2)*180/pi; 

Dss=-(Zw/Zdlt)*tan(acos(H*UA2))*  1 80/pi; 

Z=(Dss+K  1  ( 1  )*acos(H*UA2)+K  1  (2)*U*tan(acos(H*UA2)))/. 

(-Kl(4)), 

%  ESTABISHES  VEHICLE  MASS  MATRIX 

M=[l  0  0  0,0  (m-Zwd)  -Zqd  0  ;0  -Mwd  (Iy-Mqd)  0  ,0  0  0  1]; 

%  SET  TIME  LIMITS  AND  INITIAL  CONDITIONS  FOR  ODE  SOLUTION 

t0=0;tfl=2500*L/U, 
x0=[0  0  0  1], 

%  SOLVES  VEHICLE  SIMULATION  ODES 
[tl,Xl  ]=ode45('thesisde',t0,tfl  ,x0); 

%  REDEFINE  THE  OUTPUT  VECTOR  ELEMENTS  AND  FORM  DIVE  PLANE 
RESPONSE 

wl=Xl(:,2);ql=Xl(:,3);zl=Xl(:,4);thl=Xl(:,l), 

d=-(Kl(l)*thl+Kl(2)*wl+Kl(3)*ql+Kl(4)*zl); 

ll=length(d); 

forn=l:ll, 

ifd(n)>4, 

d(n)=4, 
elseifd(n)<-4, 

d(n)=-4; 
end, 
end; 
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%  PLOT  THE  RESULTS 
subplot(211); 

plot(t  1  *U/L,thl  *  1 80/pi,t  1  *U/L,d*  1 80/pi,'--');title(THET A/DELTA  vs.  TIME  '); 
xlabel('non-dimensional  time  '),ylabel('thetaydelta  (degrees)'),grid, 
gtext('~delta'),gtext(['for  a  speed  of  num2str(U) '  fps']); 

plot(tl*U/L.zl/L);title('DEPTH  vs.  TIME'),xJabel('non  dimensional  time'); 

ylabel('non  dimensional  depth  '),grid; 

pause, 

subplot(lll); 

%  THIS  IS  PROGRAM  THESISDE  M  IT  CONTAINS  THE  SUBOFF  MODEL 

SYSTEM  OF 

%  EQUATIONS  TO  BE  SOLVED  BY  PROGRAM  TFIESIS1  M 

%  SET  SYSTEM  OF  EQUATIONS  FUNCTION 

function  Xldot=thesisde(tl,Xl), 

%  SET  CONTROL  LAW 

dl=-Kl(l)*Xl(l)-Kl(2)*Xl(2)-Kl(3)*Xl(3)-Kl(4)*Xl(4); 

%  APPLY  DIVE  PLANE  SATURATION 

ifdl>4, 

dl=4; 
elseif  dl<-  4, 

dl=-4; 
end, 

%  INCLUDE  DRAG  TERMS  AND  CALCULATE  AREA  AND  CENTROID 
CD=0; 

ifXl(2)==0, 
Xl(2)-le-5, 
end, 

ifXl(3)==0, 
Xl(3)=le-5; 
end, 

Zval=bm  *((Xl(2)-xL.*Xl(3))  A3)./(abs(Xl(2)-xL.*Xl(3))), 
MvaNbm.*(((Xl(2)-xL.*Xl(3)).A3).*xL)./(abs(Xl(2)-xL.*Xl(3))); 

ZDragval=0,MDragval=0; 
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%  trapaziodal  integration 
for  n=l:length(xL)-l, 

Zdragval=0. 5  *(Zval(n)+Zval(n+ 1  ))*(xL(n+ 1  )-xL(n)); 
Mdragval=0. 5  *(Mval(n)+Mval(n+ 1  ))*(xL(n+ 1  )-xL(n)); 
ZDragval=ZDragval+Zdragval; 
MDragval =MDragval +Mdragval , 

end, 

Zdr=(0.5)*rho*CD*ZDragval; 
Mdr=(0.5)*rho*CD*MDragval; 

%  VEHICLE  SYSTEM  OF  EQUATIONS 

F1=[X1(3);  m*(U*Xl(3)+zg*Xl(3)A2)+Zq*U*Xl(3)+Zw*U*Xl(2)+UA2*Zdlt*dl-Zdr; 

m*(-zg*Xl(2)*Xl(3))+Mq*U*Xl(3)+Mw*U*Xl(2)- 
B 1  *zgb*  sin(X  1  ( 1  ))+UA2*Mdlt*d  1  -Mdr, 

-U*sin(Xl(l))+Xl(2)*cos(Xl(l))]; 

mpr=inv(M); 
Xldot=mpr*Fl; 
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C      PROGRAM  STEADY 

C 

C      STEADY  STATE  SOLUTIONS  OF  SUBMARINES 

C      IN  THE  DIVE  PLANE  AT  LOW  SPEEDS 

C 

C      INPUT  DATA 

C 

C      ICON  =  1    COMPUTATION  OF  SOLUTION  SETS  (S.S) 

C  2    COMPUTATION  OF  BIFURCATION  GRAPHS  (B.G.) 

C 

C      IV AR  =  1  :  Fn  VARIATION 

C  2  :  Xgb  VARIATION 

C 

REAL  MW,LENGTH,MASS,MDS,MDB,MD,LAMBDA 
C 

DIMENSION  B(25),X(25),BX(25) 
C 

C      GEOMETRIC  PROPERTIES  AND  HYDRODYNAMIC  COEFFICIENTS 
C 

RHO    =194 

LENGTH=  13.9792 

WRITE  (*,*)  '  ENTER  CDZ* 

READ  (V)CDZ 

CDZ    =  CDZ*0  5*RHO 

BUO=  1556  2363 

WRITE(V)  'ENTER  VALUE  OF  EXCESS  BUO' 

READ(*,*)FAC 

WEIGHT=FAC*BUO 

MASS  =  WEIGHT/32.2 

XB     =0.0 

ZB     =0.0 

WRITE  (*,*)  'ENTER  THE  METACENTRIC  HEIGHT  ZGB' 

READ(*,*)ZG 

ZW     =-0  013910*0  5*RHO*LENGTH**2 

ZDS    =-0.005603*0  5*RHO*LENGTH**2 

ZDB    =-0  005603 *0.25*RHO*LENGTH**2 

MW     =  0  010324*0.5*RHO*LENGTH**3 

MDS    =-0.002409*0  5*RHO*LENGTH**3 

MDB    =  0  002409*0  25*RHO*LENGTH**3 
C 

WRITE  (*,*)  'ENTER  THE  VALUE  OF  ALPHA' 

READ  (*,*)  ALPHA 

ZD=ZDS+ALPHA*ZDB 
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MD=MDS+ALPHA*MDB 

C 

C      DEFINE  THE  LENGTH  X  AND  BREADTH  B  TERMS  FOR  THE 

INTEGRATION 

C 

X(  1)=20.4167 
X(  2)=20.4000 
X(  3)=20.3 
X(  4)=20.2 
X(5)=20.1 
X(  6)=20.0 
X(7)=19.0 
X(8)=18.0 
X(9)=17.0 
X(10)=160 
X(  1 1 )= 1 5  1429 
X(12)=10.0 
X(13)=7.7143 
X(14)=4.0 
X(15)=3.0 
X(16)=2.0 
X(17)=l  0 
X(18)=0.7 
X(19)=06      . 
X(20)=0.5 
X(21)=0.4 
X(22)=0  3 
X(23)=0.2 
X(24)=0.1 
X(25)=0.0 
C 

B(  1)^0.00000 
B(2)=0  03178 
B(  3)=0.07920 
B(4)=0. 10074 
B(5)=0. 11243 
B(6)=0. 11724 
B(  7)=0.26835 
B(  8)=0.55025 
B(9)=0  81910 
B(10)=0.97598 
B(ll)=  1.00000 
B(  12)=  1.00000 
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B(13)=l. 00000 

B(14)=0.99282 
B(15)=0.94066 
B(16)=0. 84713 
B(17)=0. 70744 
B(18)=0.63514 
B(19)=0  60352 
B(20)=0.56627 
B(21)=0.52147 
B(22)=0  46600 
B(23)=039396 
B(24)=0. 29058 
B(25)=0  00000 

C 

DO  1  1=1,25 

B(I)=B(I)*  1.6667 
X(I)=(-X(I)+10.0)*LENGTH/20.0 
BX(I)=B(I)*X(I) 
1  CONTINUE 

CALL  TRAP(18,B,X,AREA) 
CALL  TRAP(18,BX,X,CH) 
XA=CH/AREA 

C      WRITE(*,*) 'AREA  EQUALS', AREA 

C      WRITE(*,*) 'XA  EQUALS  *,XA 
IFC=0 

C 

C      OPEN  RESULTS  FILES 

OPEN  (10,FILE='R10  RES',STATUS=fNEW) 
OPEN  ( 1 1  ,FILE='R1 1  RES',STATUS=fNEW) 
OPEN  (12,FILE='R12  RES',STATUS=rNEW) 
OPEN  (13,FILE='R13  RES',STATUS=fNEW) 
OPEN  (14,FILE='R14  RES',STATUS=fNEW) 
OPEN  (lS.FILE^RlS  RES',STATUS=fNEW) 
OPEN(16,FILE='R16.RES',STATUS=rNEW') 
OPEN  (20,FILE=,C20  RES',  ST ATUS=rNEW) 
OPEN(21,FILE='S21.RES',STATUS=,NEW') 
OPEN(22.FILE=,S22.RES,,STATUS=,NEW,) 
OPEN  (23,FILE='S23  RES*,STATUS=fNEW) 
OPEN(24,FILE=*S24.RES',STATUS=rNEW') 
OPEN  (3 1,FILE='S3 1  RES,,STATUS=rNEW) 
OPEN  (32,FILE='S32  RES,,STATUS=fNEW') 
OPEN  (33,FILE=,S33  RES",STATUS=rNEW') 
OPEN  (34,FILE='S34  RES', ST ATUS=rNEW) 
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OPEN  (41,FILE='C41  .RES,,STATUS='NEW,) 
OPEN  (42,FILE='C42  RES', ST ATUS='NEW) 
OPEN  (43,FILE='C43  RES,,STATUS=rNEW) 

C 

C      READ  DATA  INTERACTIVELY 

C 

WRITE  (*,  1001) 

READ  (*,*)     ICON 

IF  (ICON  EQ  2)  GO  TO  331 

WRITE  (*,  1009) 

READ  (*,*)     IVAR 

GO  TO  (33 1,332),  IVAR 

331  WRITE  (*,  1002) 
READ  (*,*)     UMIN 
WRITE  (*,  1003) 
READ  (*,*)     UMAX 
WRITE  (*,  1004) 
READ  (*,*)     ITER 
JU=ITER 

IF  (ICON  EQ  2)  GO  TO  332 
WRITE  (*,  1010) 
READ  (*,*)     XG 
XG=XG*  LENGTH 
GO  TO  300 

332  WRITE  (*,  1005) 
READ  (*,*)     XGMTN 
WRITE  (*,  1006) 
READ  (*,*)     XGMAX 
WRITE  (*,  1004) 
READ  (*,*)     ITER 
JXG=ITER 

XGMTN=XGMIN*LENGTH 
XGMAX=XGMAX*LENGTH 
IF  (ICON  EQ  2)  GO  TO  334 
WRITE  (*,  1012) 

READ  (*,*)     U 
334  CONTINUE 

GO  TO  300 
300  WRITE  (*,  1013) 

READ  (*,*)     ISET 

IF  (ICON  EQ  2)  GO  TO  370 

GO  TO  (350,360,400),  ISET 
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c 

C      EXACT  SOLUTION  SET 
C 

350  DO  2  1=1, ITER 

WRITE  (*,2000)  I,ITER 
GO  TO  (351,352),  IV AR 

351  U=UMTN+(UMAX-UMTN)*(I-1)/(ITER-1) 
GO  TO  355 

352  XG=XGMTN+(XGMAX-XGMIN)*(I- 1  )/(ITER- 1 ) 
355    A1=ZW*MD-MW*ZD 

A2=-CDZ*AREA*(MD-XA*ZD) 
A3=-(XG-XB)*WEIGHT*ZD 
A4=  (ZG-ZB)*WEIGHT*ZD 
CALL 
EXS0LS(IVARXL,IFC,U,XG,A1,A2,A3,A4,CDZ,AREA,ZW,ZD,ZG,MW, 
&  MD,BUO,XA) 

IF(LLGT  1)  EFC=1 

2  CONTINUE 
GO  TO  5000 

C 

C      PITCHFORK  SOLUTION  SET 

C 

360  DO  3  1=1, ITER 

WRITE  (*, 2000)  LITER 
GO  TO  (361,362),  IV AR 

361  U=UMIN+(UMAX-UMTN)*(I-1)/(ITER-1) 
GO  TO  365 

362  XG=XGMIN+(XGMAX-XGMTN)*(I- 1  )/(ITER- 1 ) 
365    A1=ZD*(ZG-ZB)*  WEIGHT 

A2=  2  0*U*U*U*CDZ*AREA*(MD-XA*ZD) 
A3=-2*U*U*(U*U*(ZW*MD-MW*ZD)+ZD*(ZG-ZB)*WEIGHT) 
A4=  2.0*U*U*ZD*(XG-XB)*WEIGHT 
A5=-ZD*(XG-XB)*WEIGHT 

CALLPITSLS(IVAR,LL,IFC,U,XG,A1,A2,A3,A4,A5) 
IF(LL.GT  1)IFC=1 

3  CONTINUE 
GO  TO  5000 

C 

C      SIMPLE  PITCHFORK  SOLUTION  SET 
C 
400  DO  401  1=1, ITER 

WRITE  (*,2000)  LITER 

GO  TO  (402,403),  IVAR 


100 


402  U=UMIN+(UMAX-UMIN)*(I-1)/(ITER-1) 
GO  TO  405 

403  XG=XGMIN+(XGMAX-XGMIN)*(I- 1  )/(ITER- 1 ) 
405    A1=ZD*(ZG-ZB)*WEIGHT 

A3=-2*U*U*(U*U*(ZW*MD-MW*ZD)+ZD*(ZG-ZB)*WEIGHT) 
A4=  2.0*U*U*ZD*(XG-XB)*WEIGHT 
CALL  SIMSLS(IVAR,LL,IFC,U,XG,A1,A3,A4) 
IF(LL.GT  1)  EFC=1 
401  CONTINUE 

GO  TO  5000 
370  GO  TO  (380,390,410,430),  ISET 
C 

C      EXACT  BIFURCATION  SET 
C 

380  DO  4  1=1, JXG 

WRITE  (*,2000)  I,JXG 

XG=XGMIN+(XGMAX-XGMIN)  *  (I- 1 )/( JXG- 1 ) 

IHS=0 

DO  5  J=1,JU 

U=UMIN+(UMAX-UMIN)*(J- 1 )/( JU- 1 ) 

A1=ZW*MD-MW*ZD 

A2=-CDZ*AREA*(MD-XA*ZD) 

A3=-(XG*WEIGHT-XB*BUO)*ZD+(WEIGHT-BUO)*MD 

A4=  (ZG*WEIGHT-ZB*BUO)*ZD 

CALL  EXNUMS(NUM,A1,A2,A3,A4,U) 

IF(JNE.1)GOT0  381 

UO=U 

NUMO=NUM 

GO  TO  5 

381  UN=U 
NUMN=NUM 

NDIF=I  AB  S(NUMO-NUMN) 

IF  (NDIF.EQ.2)  GO  TO  382 

UO=UN 

NUMO=NUMN 

GO  TO  5 

382  Ul=UO 
U2=UN 

DO  61  JJ=1,JU 
U=U1+(U2-U1)*(JJ-1)/(JU-1) 
A1=ZW*MD-MW*ZD 
A2=-CDZ*AREA*(MD-XA*ZD) 
A3=-(XG*WEIGHT-XB*BUO)*ZD+(WEIGHT-BUO)*MD 
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A4=  (ZG*WEIGHT-ZB*BUO)*ZD 

CALL  EXNUMS(NUM,A1,A2,A3,A4,U) 

IF(JJNE  1)  GO  TO  383 

UO=U 

NUMO=NUM 

GO  TO  61 

383  UN=U 
NUMN=NUM 

NDIF=IABS(NUMO-NUMN) 
IF  (NDIF.EQ.2)  GO  TO  384 
UO=UN 
NUMO=NUMN 

61  CONTINUE 

GO  TO  5 

384  Ul=UO 
U2=UN 

DO  62  JJ=1,JU 
U=U1+(U2-U1)*(JJ-1)/(JU-1) 
A1=ZW*MD-MW*ZD 
A2=-CDZ*AREA*(MD-XA*ZD) 

A3=-(XG*WEIGHT-XB*BUO)*ZD+(WEIGHT-BUO)*MD 
A4=  (ZG*WEIGHT-ZB*BUO)*ZD 
CALL  EXNUMS(NUM,A1,A2,A3,A4,U) 
IF(JJNE.l)GOT0  386 
UO=U 

NUMO=NUM 
GO  TO  62 
386        UN=U 

NUMN=NUM 

NDIF=I  AB  S(NUMO-NUMN) 

IF  (NDIF.EQ.2)  GO  TO  385 

UO=UN 

NUMO=NUMN 

62  CONTINUE 

GO  TO  5 

385  UP=(UO+UN)/2.0 
IHS=IHS+1 

UPL=SQRT(UP/(32.2*LENGTH)) 
XGL=XG/LENGTH 
IPF=10+IHS 

WRITE  (IPF,100)  XGL,UP 

UO=UN 

NUMO=NUMN 
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5    CONTINUE 

4  CONTINUE 
GO  TO  5000 
C 

C      PITCHFORK  BIFURCATION  SET 
C 

390  DO  6  1=1, JXG 

WRITE  (*,2000)  LJXG 

XG=XGMIN+(XGMAX-XGMIN)*(I- 1  )/(JXG- 1 ) 

IHS=0 

D0  7J=1,JU 

U=UMIN+(UMAX-UMIN)*(J- 1  )/(JU- 1 ) 

Al=  ZD*(ZG-ZB)*  WEIGHT 

A2=  2  0*U*U*U*CDZ*AREA*(MD-XA*ZD) 

A3=-2*U*U*(U*U*(ZW*MD-MW*ZD)+ZD*(ZG-ZB)*WEIGHT) 

A4=  2  0*U*U*ZD*(XG-XB)*WEIGHT 

A5=-ZD*(XG-XB)*WEIGHT 

CALL  PINUMS(NUM,A1,A2,A3,A4,A5,U) 

IF  (J  NE  1)  GO  TO  391 

UO=U 

NUMO=NUM 

GO  TO  7 

391  UN=U 
NUMN=NUM 

NDIF=IABS(NUMO-NUMN) 
IF  (NDIF.EQ.2)  GO  TO  392 
UO=UN 
NUMO=NUMN 

GO  TO  7 

392  Ul=UO 
U2=UN 

DO  71  JJ=1,JU 
U=U1+(U2-U1)*(JJ-1)/(JU-1) 
Al=  ZD*(ZG-ZB)*  WEIGHT 
A2=  2  0*U*U*U*CDZ*AREA*(MD-XA*ZD) 
A3=-2*U*U*(U*U*(ZW*MD-MW*ZD)+ZD*(ZG-ZB)*WEIGHT) 
A4=  2  0*U*U*ZD*(XG-XB)*WEIGHT 
A5=-ZD*(XG-XB)*WEIGHT 
CALL  PINUMS(NUM,ALA2,A3,A4,A5,U) 
IF(JJNE.l)GOT0  393 
UO=U 

NUMO=NUM 
GO  TO  71 
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393  UN=U 
NUMN=NUM 

NDIF=I  AB  S(NUMO-NUMN) 
IF  (NDIF.EQ.2)  GO  TO  394 
UO=UN 
NUMO=NUMN 

71  CONTINUE 

STOP  1111 

394  Ul=UO 
U2=UN 

D0  72JJ=1JU 
U=U1+(U2-U1)*(JJ-1)/(JU-1) 
Al=  ZD*(ZG-ZB)*  WEIGHT 
A2=  2  0*U*U*U*CDZ*AREA*(MD-XA*ZD) 
A3=-2*U*U*(U*U*(ZW*MD-MW*ZD)+ZD*(ZG-ZB)*WEIGHT) 
A4=  2  0*U*U*ZD*(XG-XB)*WEIGHT 
A5=-ZD*(XG-XB)*WEIGHT 
CALL  PINUMS(NUM,A1,A2,A3,A4,A5,U) 
IF(JJNE  1)  GO  TO  396 
UO=U 

NUMO=NUM 
GO  TO  72 
396        UN=U 

NUMN=NUM 

NDIF=IABS(NUMO-NUMN) 

IF  (NDIF.EQ.2)  GO  TO  395 

UO=UN 

NUMO=NUMN 

72  CONTINUE 

STOP  1112 

395  UP=(UO+UN)/2.0 
IHS=IHS+1 

UPL=SQRT(UP/(32.2*LENGTH)) 
XGL=XG/LENGTH 
IPF=10+IHS 

WRITE  (IPF,100)  XGL,UP 
UO=UN 
NUMO=NUMN 
7    CONTINUE 
6  CONTINUE 
GO  TO  5000 
C 
C      SIMPLE  PITCHFORK  BIFURCATION  SET 
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c 

410  DO  411  I=1,JXG 

WRITE  (*,2000)  I,JXG 

XG=XGMIN+(XGMAX-XGMIN)*(I- 1 )/( JXG- 1 ) 

IHS=0 

DO  412  J=1,JU 

U=UMIN+(UMAX-UMIN)*(J- 1  )/(JU- 1 ) 

Al=  ZD*(ZG-ZB)*WEIGHT 

A3=-2*U*U*(U*U*(ZW*MD-MW*ZD)+ZD*(ZG-ZB)*WEIGHT) 

A4=  2  0*U*U*ZD*(XG-XB)*WEIGHT 

CALL  SINUMS(NUM,A1,A3,A4,U) 

IF(JNE.l)GOT0  413 

UO=U 

NUMO=NUM 

GO  TO  412 

413  UN=U 
NUMN=NUM 

NDIF=IABS(NUMO-NUMN) 
IF  (NDIF.EQ.2)  GO  TO  414 
UO=UN 
NUMO=NUMN 

GO  TO  412 

414  Ul=UO 
U2=UN   • 

DO  415  JJ=1,JU 
U=U1+(U2-U1)*(JJ-1)/(JU-1) 
Al=  ZD*(ZG-ZB)*WEIGHT 

A3=-2*U*U*(U*U*(ZW*MD-MW*ZD)+ZD*(ZG-ZB)*WEIGHT) 
A4=  2.0*U*U*ZD*(XG-XB)*WEIGHT 
CALL  SINUMS(NUM,A1,A3,A4,U) 
IF(JJNE.l)GOT0  416 
UO=U 

NUMO=NUM 
GOT0415 

416  UN=U 
NUMN=NUM 

NDIF=IABS(NUMO-NUMN) 
IF  (NDIF.EQ.2)  GO  TO  417 
UO=UN 
NUMO=NUMN 

415  -  CONTINUE 

STOP  1111 

417  Ul=UO 
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U2=UN 

DO  418  JJ=1,JU 

U=U1+(U2-U1)*(JJ-1)/(JU-1) 

Al=  ZD*(ZG-ZB)*  WEIGHT 

A3=-2*U*U*(U*U*(ZW*MD-MW*ZD)+ZD*(ZG-ZB)*WEIGHT) 

A4=  2  0*U*U*ZD*(XG-XB)*WEIGHT 

CALL  SINUMS(NUM,A1,A3,A4,U) 

IF(JJ.NE.l)GOT0  419 

UO=U 

NUMO=NUM 

GO  TO  418 

419  UN=U 
NUMN=NUM 

NDIF=I  AB  S(NUMO-NUMN) 
IF  (NDIF.EQ.2)  GO  TO  420 
UO=UN 
NUMO=NUMN 
418      CONTINUE 
STOP  1112 

420  UP=(UO+UN)/2.0 

IHS=IHS+1 

UPL=SQRT(UP/(32.2*LENGTH)) 
XGL=XG/LENGTH 
IPF=10+IHS 

WRITE  (IPF,100)  XGL,UP 
UO=UN 
NUMO=NUMN 
412    CONTINUE 
411  CONTINUE 
GO  TO  5000 
C 

C      ANALYTIC  BIFURCATION  SET 
C 
430  WRITE  (*,  101 7) 
READ  (*,*)     ICUSP 
GO  TO  (432,433),  ICUSP 
C 

C         EXACT"  CUSP 
C 
432  DO  431  I=1,JU 

-  WRITE  (*, 2000)  I, JU 
U=UMIN+(UMAX-UMIN)*(I- 1  )/(JU- 1 ) 
DELTA=32  2*(ZW*MD-MW*ZD)/(ZD*WEIGHT) 
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LAMBDA=U*U/(32  2*(ZG-ZB)) 

BETA2=(8.0*U*U*(1  0+DELTA*LAMBDA)**3)/(27  0*LAMBDA**2) 

IF  (BETA2.LT.0.0)  GO  TO  43 1 

BETA=SQRT(BETA2) 

XGB=U*U*BETA/32  2 

WRITE  (10,100)  U,XGB/LENGTH 
43 1  CONTINUE 
GO  TO  5000 
C 

C        APPROXIMATE"  CUSP 
C 

433  DO  434  1=1, JU 

WRITE  (*,2000)  I,JU 

U=UMTN+(UMAX-UMTN)*(I- 1  )/(JU- 1 ) 

DELTA=32  2*(ZW*MD-MW*ZD)/(ZD*WEIGHT) 

LAMBDA=U*U/(32.2*(ZG-ZB)) 

BETA2=(8.0*U*U*(1  0+DELTA*LAMBDA)**3)/27.0 

IF  (BETA2  LT.0.0)  GO  TO  434 

BETA=SQRT(BETA2) 

XGB=U*U*BETA/32  2 

WRITE  (10,100)  U,XGB/LENGTH 

434  CONTINUE 
C 

5000  STOP 

1001  FORMAT  ('  ENTER  1  :  SOLUTION  SETS        ',/ 
1  2  :  BIFURCATION  GRAPHS    ') 

1009  FORMAT  ('  ENTER  1  :  U  VARIATION  ',/ 
1  '        2    XG  VARIATION ') 

1002  FORMAT  ('  ENTER  MTNTMUM  VALUE  OF  U) 

1003  FORMAT  ('  ENTER  MAXIMUM  VALUE  OF  U) 

1004  FORMAT  ('  ENTER  NUMBER  OF  INCREMENTS') 

1005  FORMAT  ('  ENTER  MINIMUM  VALUE  OF  XG/L') 

1006  FORMAT  ('  ENTER  MAXIMUM  VALUE  OF  XG/L') 

1010  FORMAT  ('  ENTER  VALUE  OF  XG/L') 

1012  FORMAT  ('  ENTER  VALUE  OF  U) 

1013  FORMAT  ('  ENTER  1  :  EXACT  COMPUTATION       ',/ 


1  '        2 

2  '        3 

3  '        4 


PITCHFORK  APPROXIMATION',/ 
SIMPLE  PITCHFORK        *,/ 
ANALYTIC  SET  ') 

1017  FORMAT  ('  ENTER  1  :  EXACT  CUSP  ',/ 
1  2  :  APPROXIMATE  CUSP  ') 

2000  FORMAT  (215) 

100  FORMAT  (2E 15  5) 
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END 

C 

SUBROUTINE  TRAP(N,A,B,OUT) 
C 

C      NUMERICAL  INTEGRATION  ROUTINE  USING  THE  TRAPEZOIDAL  RULE 
C 

DIMENSION  A(1),B(1) 

N1=N-1 

OUT=0.0 

DO  1  I=1,N1 

OUT  1=0.5  *(  A(I)+A(I+ 1  ))*(B(I+ 1  )-B(I)) 
OUTOUT+OUT1 
1  CONTINUE 

RETURN 

ENT) 
C - 

SUBROUTINE  EXSOLS(IVAR,L,IFC,U,XG,Al,A2, A3, A4,CDZ, ARE A,ZW, 

&ZD,ZG,MW,MD,BUO,XA) 
C 

C      IT  COMPUTES  AND  PRINTS  THE  EXACT  SOLUTION  SET  THETA 
C      (DEG)  VERSUS  U  OR  XG 
C 

REAL  MW,MD,PRNT,PI 

DIMENSION  VF(5,2) 

PI=4.0*ATAN(1.0) 

IFdVAR.EQ  1)PRNT=U 

DF  (IVAR.EQ.2)  PRNT=XG 
C 

C      FIND  FIRST  ESTIMATE  OF  SOLUTIONS 
C 

L=0 

VA=-89.5 

VA=VA*PI/180  0 

VAO=VA 

V0=THETEQ(1,VA,A1,A2,A3,A4,U) 

VAMIN=-89.5 

VAMAX=+89.5 

IVA=5000 

DO  10I=2,IVA 
AI=I 

VA=VAMIN+(VAMAX-VAMIN)*(I-1)/(IVA-1) 
C         VA=AI-90  5 

VA=VA*Py  180.0 
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VAN=VA 

VN=THETEQ(1,VA,A1,A2,A3,A4,U) 
VP=VO*VN 

IF  (VP  GEO  0)  GO  TO  11 
L=L+1 

VF(L,l)=VAO 
VF(L,2)=VAN 
11    VO=VN 

VAO=VAN 
10  CONTINUE 
C 

C      EXACT  COMPUTATION  OF  SOLUTIONS  VIA  NEWTON'S  METHOD 
C 

E=l.E-8 

IEND=500 

DO20J=l,L 

X=(VF(J,l)+VF(J,2))/2  0 
C 

X1=X 
IXX=0 

IF(IXXEQ0)GOTO35 
C 

F=THETEQ(  1  ,X,  Al ,  A2,  A3,  A4,U) 
FDER=THETEQ(2,X,A1,A2,A3,A4,U) 
DO30K=l,IEND 
IF  (FDER.EQ.O  0)  STOP  1001 
DX=F/FDER 
X1=X-DX 

F=THETEQ(  1  ,X1 ,  Al ,  A2,  A3,  A4,U) 
FDER=THETEQ(2,X1,A1,A2,A3,A4,U) 
IF  (FEQ  0  0)  GO  TO  35 
A=ABS(X1-X) 
IF  (A-E)  35,35,40 
40      X=X1 
30    CONTINUE 
GO  TO  20 
35    SOL=X1*180/PI 
C        WRITE  (*,*)  XI 
JPR=10+J 
WD=U*TAN(X1) 

DLT=((CDZ*AREA*WD*ABS(WD))-(ZW*U*WD))/(ZD*U*U) 
IF  ((ABS(DLT)  GT  0  4)  AND  (CDZ  NE.0.0))  THEN 
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CALL 
SATUR(U,XG,CDZ,AREA,ZW,ZD,ZG,MW,MD,BUO,XA,DLT,PRNT) 
ELSEIF  ((ABS(DLT).GT  0  4).AND  (CDZ  EQ  0  0))  THEN 
CALL 
SATUR1(U,XG,CDZ,AREA,ZW,ZD,ZG,MW,MD,BU0,XA,DLT,PRNT) 
ENDIF 
IF  (L.EQ  1  AND  IFC  EQ  0)  WRITE  (10,100)  PRNT,SOL,DLT*180/PI 
IF  (L  EQ  1  AND  IFC  EQ  1 )  WRITE  ( 1 6, 1 00)  PRNT,SOL,DLT*  1 80/PI 
IF(LGT  1)  WRITE  ( JPR,1 00)  PRNT,SOL,DLT*l  80/PI 

20  CONTINUE 

RETURN 
100  FORMAT  (3E1 5. 5) 
END 

C - 

SUBROUTINE 
SATUR(U,XG,CDZ,AREA,ZW,ZD,ZG,MW,MD,BUO,XA,DLT,PRNT) 

REAL  MW,MD,PRNT 
C 

C      USED  TO  COMPUTE  SATURATION  CONDITIONS  FOR  THE  DIVE  PLANE 
WITH  DRAG 

C      COEFFICENT  INCLUDED 
C 

CD=CDZ 

IF  (DLT  GT  0  0)  THEN 

Wl=(-ZW*U+SQRT((ZW*U)*(ZW*U)-4*(-CD*AREA)*ZD*U*U*0.4)) 
&  /(2*(-CD*AREA)) 

IF(W1.GT0  0)THEN 
A=MW*U*W1-CD*XA*AREA*W1*W1-XG*BUO+MD*U*U*0  4 
B=-20*ZG*BUO 

C=MW*U*W1-CD*XA*AREA*W1*W1+XG*BUO+MD*U*U*0.4 
X 1  =(-B+SQRT((B  1  *B  1  )-4*(A*C)))/(2.0*  A) 
SOL1=2.0*ATAN(X1) 
X2=(-B-SQRT((B  1  *B  1  )-4*(  A*C)))/(2.0*  A) 
SOL2=2.0*ATAN(X2) 
WRITE  (21,200)  PRNT,SOLl,SOL2 
ENDIF 
W2=(-ZW*U-SQRT((ZW*U)*(ZW*U)-4*(-CD*AREA)*ZD*U*U*0.4)) 
&  /(2*(-CD*AREA)) 

IF  (W2  GT.0.0)  THEN 
A=MW*U*W2-CD*XA*AREA*W2*W2-XG*BUO+MD*U*U*0  4 

B=-20*ZG*BUO 
C=MW*U*W2-CD*XA*AREA*W2*W2+XG*BUO+MD*U*U*0  4 

X1=(-B+SQRT((B  1  *B  1  )-4*(A*C)))/(2.0*  A) 
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SOL1=2.0*ATAN(X1) 

X2=(-B-SQRT((B  1  *B  1  )-4*(  A*C)))/(2.0*  A) 

SOL2=2.0*ATAN(X2) 

WRITE  (22,200)  PRNT,SOLl,SOL2 
ENDIF 
W3=(-ZW*U+SQRT((ZW*U)*(ZW*U)-4*(CD*AREA)*ZD*U*U*0  4)) 
&  /(2*(CD*AREA)) 

IF(W3.LT.0.0)THEN 

A=MW*U*W3+CD*XA*AREA*W3*W3-XG*BUO+MD*U*U*0  4 

B=-20*ZG*BUO 

C=MW*U*W3+CD*XA*AREA*W3*W3+XG*BUO+MD*U*U*0  4 

X  1=(-B+SQRT((B  1  *B  1  )-4*(  A*C)))/(2.0*  A) 

SOL  1=2.0*  AT  AN(X1) 

X2=(-B-SQRT((B  1  *B  1  )-4*(A*C)))/(2.0*  A) 

SOL2=2.0*ATAN(X2) 

WRITE  (23,200)  PRNT, SOLI, SCO 
ENDIF 
W4=(-ZW*U-SQRT((ZW*U)*(ZW*U)-4*(CD*AREA)*ZD*U*U*0  4)) 
&  /(2*(CD*AREA)) 

IF  (W4.LT.0.0)  THEN 

A=MW*U*W4+CD*XA*AREA*W4*W4-XG*BUO+MD*U*U*0  4 

B=-20*ZG*BUO 

C=MW*U*W4+CD*XA*AREA*W4*W4+XG*BUO+MD*U*U*0.4 

X1=(-B+SQRT((B1*B1)-4*(A*C)))/(2.0*A) 

SOL1=2.0*ATAN(X1) 

X2=(-B-SQRT((B  1  *B  1  )-4*(A*C)))/(2.0*  A) 

SOL2=2.0*ATAN(X2) 

WRITE  (24,200)  PRNT, SOL  1,S0L2 
ENDIF 
ELSE 

Wl=(-ZW*U+SQRT((ZW*U)*(ZW*U)-4*(-CD*AREA)*ZD*U*U* 
&  (-0  4)))/(2*(-CD*AREA)) 

IF(W1.GT.0.0)THEN 

A=MW*U*W1-CD*XA*AREA*W1*W1-XG*BUO+MD*U*U*(-0  4) 

B=-20*ZG*BUO 

C=MW*U*W1-CD*XA*AREA*W1*W1+XG*BUO+MD*U*U*(-0  4) 

X1=(-B+SQRT((B  1  *B  1  )-4*(A*C)))/(2  0* A) 

SOL1=2.0*ATAN(X1) 

X2=(-B-SQRT((B  1  *B  1  )-4*( A*C)))/(2  0*  A) 

SOL2=2.0*ATAN(X2) 

WRITE  (31,200)  PRNT,  SOL  1,S0L2 
ENDIF 
W2=(-ZW*U-SQRT((ZW*U)*(ZW*U)-4*(-CD*AREA)*ZD*U*U* 
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&  (-0.4)))/(2*(-CD*AREA)) 

IF  (W2.LT.0.0)  THEN 

A=MW*U*W2-CD*XA*AREA*W2*W2-XG*BUO+MD*U*U*(-0  4) 

B=-20*ZG*BUO 

C=MW*U*W2-CD*XA*AREA*W2*W2+XG*BUO+MD*U*U*(-0.4) 

X 1  =(-B+SQRT((B  i  *B  1  )-4*(A*C)))/(2.0*  A) 

SOL1=2.0*ATAN(X1) 

X2=(-B-SQRT((B  1  *B  1  )-4*(A*C)))/(2.0*  A) 

SOL2=2.0*ATAN(X2) 

WRITE  (32,200)  PRNT,SOL21,SOL22 
ENDIF 
W3=(-ZW*U+SQRT((ZW*U)*(ZW*U)-4*(CD*AREA)*ZD*U*U* 
&  (-0  4)))/(2*(CD*AREA)) 

IF  (W3  LT  0  0)  THEN 

A=MW*U*W3+CD*XA*AREA*W3*W3-XG*BUO+MD*U*U*(-0.4) 

B=-20*ZG*BUO 

C=MW*U*W3+CD*XA*AREA*W3*W3+XG*BUO+MD*U*U*(- 


0.4) 

X 1  =(-B+SQRT((B  1  *B  1  )-4*(  A*C)))/(2.0*  A) 
SOL1=2  0*ATAN(X1) 
X2=(-B-SQRT((B  1  *B  1  )-4*(A*C)))/(2.0*  A) 
SOL2=2.0*ATAN(X2) 
WRITE  (33,200)  PRNT,SOLl,SOL2 
ENDIF 
W4=(-ZW*U+SQRT((ZW*U)*(ZW*U)-4*(CD*AREA)*ZD*U*U* 
&  (-0  4)))/(2*(CD*AREA)) 

IF  (W4.LT.0.0)  THEN 
A=MW*U*W4+CD*XA*AREA*W4*W4-XG*BUO+MD*U*U*(-0.4) 
B=-2.0*ZG*BUO 

C=MW*U*W4+CD*XA*AREA*W4*W4+XG*BUO+MD*U*U*(- 
0.4) 

X1=(-B+SQRT((B  1  *B  1  )-4*(A*C)))/(2.0*A) 
SOL1=2  0*ATAN(X1) 
X2=(-B-SQRT((B  1  *B  1  )-4*(  A*C)))/(2.0*  A) 
SOL2=2.0*ATAN(X2) 
WRITE  (34,200)  PRNT,SOLl,SOL2 
ENDIF 
ENDIF 
200  FORMAT  (3E15 .5) 
END 

C - 

SUBROUTINE 
SATURl(U,XG,CDZ,AREAZW,ZD,ZG,MW,MD,BUO,XA,DLT,PRNT) 
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REAL  MW,MD,PRNT 
C 

C      COMPUTES  DIVE  PLANE  SATURATION  NEGLECTING  DRAG 
C 

PI=4.0*ATAN(1.0) 

A1=MD*U*U*0  4-MW*U*U*ZD*0  4-XG*ZW*BUO 

A2=-MD*U*U*0  4+MW*U*U*ZD*0  4-XG*ZW*BUO 

B1=-2  0*ZG*ZW*BUO 

C 1  =-MW*U*U*ZD*0.4+XG*ZW*BUO+MD*U*U*ZW*0.4 

C2=MW*U*U*ZD*0  4+XG*ZW*BUO-MD*U*U*ZW*0  4 

XCDl=((-Bl/Al)+SQRT(((Bl/Al)**2)-4.0*(Cl/Al)))/2.0 

XCD2=((-Bl/Al)-SQRT(((Bl/Al)**2)-4.0*(Cl/Al)))/2.0 

XCD3=((-Bl/A2)+SQRT(((Bl/A2)**2)-4.0*(C2/A2)))/2.0 

XCD4=((-Bl/A2)-SQRT(((Bl/A2)**2)-4.0*(C2/A2)))/2.0 

XC  1=2.0*  AT  AN(XCD1)*1 80/PI 

XC2=2.0*  ATAN(XCD2)*  1 80/PI 

XC3=2  0* ATAN(XCD3)*  1 80/PI 

XC4=2.0*  ATAN(XCD4)*  1 80/PI 

WRITE  (20, 1 00)  PRNT,XC  1  ,DLT*  1 80/PI 

WRITE  (41,100)  PRNT,XC2,DLT*1 80/PI 

WRITE  (42,100)  PRNT,XC3,DLT*1 80/PI 

WRITE  (43,100)  PRNT,XC4,DLT*1 80/PI 
100  FORMAT  (3E 15. 5) 

END 
C - 

SUBROUTINE  PITSLS(IVAR,LL,IFC,U,XG,A1,A2,A3,A4,A5) 
C 

C      IT  COMPUTES  AND  PRINTS  THE  PITCHFORK  SOLUTION  SET  THETA 
C      (DEG  )  VERSUS  U  OR  XG 
C 

DIMENSION  VF(5,2) 

PI=4.0*ATAN(1.0) 

IF(IVAR.EQ  1)PRNT=U 

IF  (IVAR.EQ  2)  PRNT=XG 
C 

C      FIND  FIRST  ESTIMATE  OF  SOLUTIONS 
C 

L=0 

VA=-89  5 

VA=VA*PI/180.0 

VAO=VA 

V0=PITCEQ(1,VA,A1,A2,A3,A4,A5,U) 

DO  10  1=2,180 
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AI=I 

VA=AI-90.5 
VA=VA*PI/180.0 
VAN=VA 

VN=PITCEQ(1,VA,A1,A2,A3,A4,A5,U) 
VP=VO*VN 

IF  (VP  GEO  0)  GO  TO  11 
L=L+1 

VF(L,l)=VAO 
VF(L,2)=VAN 
1 1    VO=VN 

VAO=VAN 
10  CONTINUE 
C 

C      EXACT  COMPUTATION  OF  SOLUTIONS  VIA  NEWTON'S  METHOD 
C 

E=l.E-5 

IEND=500 

DO20J=l,L 

X=(VF(J,l)+VF(J,2))/2.0 
F=PITCEQ(1,X,A1,A2,A3,A4,A5,U) 
C 

X1=X 
IXX=0 

IF(IXXEQ0)GOTO35 
C 

FDER=PITCEQ(2,X,A1,A2,A3,A4,A5,U) 
DO30K=l,IEND 
IF  (FDER.EQ  0  0)  STOP  1001 
DX=F/FDER 
X1=X-DX 

F=PITCEQ(  1  ,X  1 ,  Al ,  A2,  A3,  A4,  A5,U) 
FDER=PITCEQ(2,X1,A1,A2,A3,A4,A5,U) 
IF(FEQ.0  0)GOTO35 
A=ABS(X1-X) 
IF  (A-E)  35,35,40 
40      X=X1 
30    CONTINUE 
GO  TO  20 
35    SOL=X1*180.0/PI 
JPR=10+J 

IF  (L.EQ.  1  AND TFC.EQ.O)  WRITE  ( 10, 100)  PRNT,SOL 
IF  (L  EQ  1  AND TFC  EQ  1 )  WRITE  (16,100)  PRNT,SOL 
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IF(L.GTl)  WRITE  ( JPR,1 00)  PRNT,  SOL 

20  CONTINUE 

RETURN 
100  FORMAT  (2E 15  5) 
END 

C 

SUBROUTINE  SIMSLS(IVAR,L,DFC,U,XG,A1,A3,A4) 
C 

C      IT  COMPUTES  AND  PRINTS  THE  SIMPLE  PITCHFORK  SOLUTION  SET 
THETA 

C      (DEG  )  VERSUS  U  OR  XG 
C 

DIMENSION  VF(5,2) 
PI=4.0*ATAN(1.0) 
IF(IVAREQ  1)PRNT=U 
IF  (IVAR.EQ.2)  PRNT=XG 
C 

C      FIND  FIRST  ESTIMATE  OF  SOLUTIONS 
C 

L=0 

VA=-89.5 
VA=VA*PI/180.0 
VAO=VA 

V0=SIMPEQ(1,VA,A1,A3,A4,U) 
DO  10  1=2,180 
AI=I 

VA=AI-90  5 
VA=VA*PI/ 180.0 
VAN=VA 

VN=SIMPEQ(1,VA,A1,A3,A4,U) 
VP=VO*VN 

IF  (VP  GEO  .0)  GO  TO  11 
L=L+1 

VF(L,l)=VAO 
VF(L,2)=VAN 
11    VO=VN 

VAO=VAN 
10  CONTINUE 
C 

C      EXACT  COMPUTATION  OF  SOLUTIONS  VIA  NEWTON'S  METHOD 
C 

E=l.E-5 
IEND=500 
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DO20J=l,L 

X=(VF(J,l)+VF(J,2))/2.0 
C 

X1=X 
IXX=0 

IF(IXXEQ0)GOTO35 
C 

F=SIMPEQ(1,X,A1,A3,A4,U) 
FDER=SIMPEQ(2,X,A1,A3,A4,U) 
DO30K=l,IEND 
IF  (FDER.EQ.O.O)  STOP  1001 
DX=F/FDER 
X1=X-DX 

F=SIMPEQ(1,X1,A1,A3,A4,U) 
FDER=SIMPEQ(2,X1,A1,A3,A4,U) 
IF  (FEQ .0.0)  GO  TO  35 
A=ABS(X1-X) 
IF  (A-E)  35,35,40 
40      X=X1 
30    CONTINUE 
GO  TO  20 
35    SOL=X1*180.0/PI 
JPR=10+J 

IF  (L  EQ  1  AND TEC  EQ  0)  WRITE  (10, 100)  PRNT,SOL 
IF(LEQ  1  ANDTPC.EQ  1 )  WRITE  ( 1 6, 1 00)  PRNT,SOL 
IF(L.GT.l)  WRITE  ( JPR,1 00)  PRNT, SOL 

20  CONTINUE 

RETURN 
100  FORMAT  (2E 15. 5) 
END 

C - - 

SUBROUTINE  EXNUMS(NUM,A1,A2,A3,A4,U) 
C 

C   IT  COMPUTES  THE  NUMBER  OF  SOLUTIONS  OF  THE  EXACT  SOLUTION 
SET 
C 

PI=4.0*ATAN(1.0) 
L=0 

VA=-89.5 
VA=VA*PI/180.0 

VO=THETEQ(  1,VA,  Al ,  A2,A3,  A4,U) 
DO  101=2,180 
AI=I 
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VA=AI-90.5 
VA=AI*PI/180.0 

VN=THETEQ(  1,VA,A1,A2,A3,A4,U) 
VP=VO*VN 

IF  (VP  GE  0.0)  GO  TO  11 
L=L+1 
11    VO=VN 

10  CONTINUE 
NUM=L 
RETURN 
END 

C - 

SUBROUTINE  PINUMS(NUM,A1,A2,A3,A4,A5,U) 
C 

C      IT  COMPUTES  THE  NUMBER  OF  SOLUTIONS  OF  THE  PITCHFORK 
SOLUTION  SET 
C 

PI=4.0*ATAN(1.0) 
L=0 

VA=-89  5 
VA=VA*PI/180.0 

V0=PITCEQ(1,VA,A1,A2,A3,A4,A5,U) 
DO  10  1=2,180 
AI=I 

VA=AI-90.5 
VA=VA*PI/180  0 

VN=PITCEQ(1,VA,A1,A2,A3,A4,A5,U) 
VP=VO*VN 

IF  (VP  GE.0.0)  GO  TO  11 
L=L+1 

1 1  VO=VN 
10  CONTINUE 

NUM=L 

RETURN 

END 
C 

SUBROUTINE  SINUMS(NUM,A1,A3,A4,U) 
C 

C      IT  COMPUTES  THE  NUMBER  OF  SOLUTIONS  OF  THE  SIMPLE  PITCHFORK 
C      SOLUTION  SET 
C 

PI=4.0*ATAN(1.0) 

L=0 
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VA=-89.5 
VA=VA*PI/180.0 
V0=SIMPEQ(1,VA,A1,A3,A4,U) 
DO  101=2,180 

AI=I 

VA=AI-90.5 

VA=VA*PI/180.0 

VN=SIMPEQ(1,VA,A1,A3,A4,U) 

VP=VO*VN 

IF  (VP  GEO  0)  GO  TO  11 

L=L+1 
1 1    VO=VN 
10  CONTINUE 
NUM=L 
RETURN 
END 

C 

FUNCTION  THETEQ(K,XA,A1,A2,A3,A4,U) 
C 

C      IT  COMPUTES  THE  VALUE  OF  THE  EXACT  EQUATION  FOR 
C      THETABAR  FOR  A  GIVEN  VALUE  OF  THETA 
C 

C      K  =  1  :  COMPUTE  THE  VALUE  OF  THE  FUNCTION 
C      K  =  2  :  COMPUTE  THE  VALUE  OF  ITS  DERIVATIVE 
C 

SP=SIN(XA) 

CP=COS(XA) 

AP=ABS(SP) 

TP=TAN(XA) 

ATP=ABS(TP) 

GO  TO  (10,20),  K 
10THETEQ=A1*U*U*SP*CP+A2*U*U*SP*AP+A3*CP**3+A4*SP*CP**2 
C     10THETEQ=A1*U*TP+A2*U*U*TP*ATP+A3*CP+A4*SP 

GO  TO  50 
20  THETEQ=Al*U*U*CP*CP-Al*U*U*SP*SP-3  0*A3*CP*CP*SP+A4*CP**3 

&       -2.0*A4*SP*SP*CP+2.0*A2*U*U*CP*AP 
50  RETURN 

END 
C 

FUNCTION  PITCEQ(K,XA,A1,A2,A3,A4,A5,U) 
C 

C      IT  COMPUTES  THE  VALUE  OF  THE  PITCHFORK  EQUATION  FOR 
C      THETABAR  FOR  A  GIVEN  VALUE  OF  THETA 
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c 

C      K  =  1  :  COMPUTE  THE  VALUE  OF  THE  FUNCTION 
C      K  =  2    COMPUTE  THE  VALUE  OF  ITS  DERIVATIVE 
C 

TP  =TAN(XA) 

CP  =COS(XA) 

CP2=CP*CP 

AP  =ABS(TP) 

GO  TO  (10,20),  K 
1 0  PITCEQ=A1  *(U*TP)*  *3+A2*U*U*TP*  AP+A3 *U*TP+A4+A5 *U*U*TP*TP 

GO  TO  50 
20PITCEQ=3.0*U*U*U*TP/CP2+2.0*A2*U*U*AP/CP2+A3*U/CP2 

&       +2.0*A5*U*U*TP/CP2 
50  RETURN 

END 
C - - 

FUNCTION  SIMPEQ(K,XA,A1,A3,A4,U) 
C 

C      IT  COMPUTES  THE  VALUE  OF  THE  SIMPLE  PITCFIFORK  EQUATION  FOR 
C      TFIETABAR  FOR  A  GIVEN  VALUE  OF  THETA 
C 

C      K  -  1    COMPUTE  THE  VALUE  OF  THE  FUNCTION 
C      K  =  2    COMPUTE  THE  VALUE  OF  ITS  DERIVATIVE 
C 

TP  =TAN(XA) 

CP  =COS(XA) 

CP2=CP*CP 

AP  =ABS(TP) 

GO  TO  (10,20),  K 
10  SIMPEQ=A1*(U*TP)**3+A3*U*TP+A4 

GO  TO  50 
20  SIMPEQ=3  0*U*U*U*TP/CP2+A3*U/CP2 
50  RETURN 

END 
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%  THIS  IS  PROGRAM  BIFSUM  IT  COMPUTES  THE  BIASED  BIFURCATION 

CUSP  FOR 

%  DELTA  SATURATED  INCLUDING  DRAG  TERMS  IT  USES  THE  FZERO 

FUNCTION  TO 

%  SOLVE  THE  EQUATIONS 

global  W  Bl  U  n  Zw  Zdlt  A  XA  CD  rho, 

%  ESTABLISH  GEOMETRIC  PARAMETERS 

W=l 556  2363 ;B1  =  1. 0001  *W; 

A=19.8473;XA=0. 20126; 

CD=0.3; 

Iy=561.32; 

g-32.2; 

m=W/g; 

rho=1.94; 

L=13.9792; 

xb=0, 

Alphas=l, 
Alphab=0; 
zg=0.1; 

zb=0; 
zgb=zg-zb; 

%  NON-DIMENSIONING  FACTORS 

ndl=.5*rho*LA2, 

nd2=5*rho*LA3; 

nd3=.5*rho*LA4, 

nd4=.5*rho*LA5, 

%  HYDRODYNAMIC  COEFFICIENTS 

Zqdnd=-6. 33e-4;Zwdnd=- 1  4529e-2;Zqnd=7. 545e-3  ;Zwnd=- 1 .39 1  e-2; 
Zds=-5.603e-3;Zdb=0.5*(-5.603e-3);Zdltnd=(Alphas*Zds+Alphab*Zdb); 
Mqdnd=-8.8e-4;Mwdnd=-5.61e-4;Mqnd=-3.702e-3;Mwnd=1.0324e-2; 
Mds=-0.002409;Mdb=0.5*(0.002409);Mdltnd=(AIphas*Mds+AJphab*Mdb); 

Zqd=nd3  *Zqdnd;Zwd=nd2*Zwdnd,Zq=nd2*Zqnd,Zw=nd  1  *Zwnd, 
Zdlt=ndl*Zdltnd, 

Mqd=nd4*Mqdnd,Mwd=nd3*Mwdnd,Mq=nd3*Mqnd,Mw=nd2*Mwnd; 
Mdlt=nd2*Mdltnd: 
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%  ESTABLISH  SPEED  RANGE 

U=9:0.01:2.6; 
forn=l:length(U), 

%  SOLVES  FORCE  EQUATION  USING  FZERO  FUNCTION 

th(n)=fzero('satur',(-Zdlt/Zw),0.00001); 

%  COMPUTES  RELATIONSHIP  BETWEEN  U  AND  Xgb  WITH  RESULTS  FROM 
ABOVE  EQUATION 

xg(n)=(Mw*U(n)A2*tan(th(n))-zg*W*sin(th(n))-(0.5)*rho*CD*A*XA*U(n)A2*.. 

tan(th(n))*abs(tan(th(n)))+Mdlt*U(n)A2*(0.4))/(W*cos(th(n))); 

end, 

XgbL=xg/l  3.9792; 

Lp=sqrt(UA2/(3.22)), 

for  n=l:length(U); 

th  1  (n)=fzero('satur  1  \(-Zdlt/Zw),0.0000 1 ); 

xgl(n)=(Mw*U(n)A2*tan(thl(n))-zg*W*sin(thl(n))-(0  5)*rho*CD*A* 

XA*U(n)A2*tan(thl(n))*abs(tan(thl(n)))+Mdlt*U(n)A2*(-0.4))/(W*cos(thl(n))), 

end, 

XgbLl=xgl/13.9792, 


%  THIS  IS  PROGRAM  SATUR.M.  IT  ESTABLISHES  THE  FUNCTION  TO  BE 

SOLVED  IN 

%  THE  BIFSU  M  PROGRAM 

function  y=satur(th); 

y=Zw*U(n)A2*tan(th)+(  W-B 1  )*cos(th)+Zdlt*U(n)A2*(0.4). . . 

-(0.5)*rho*CD*A*U(n)A2*tan(th)*abs(tan(th)); 

%  THIS  IS  PROGRAM  SATUR1.M  IT  ESTABLISHES  THE  FUNCTION  TO  BE 

SOLVED  IN 

%  THE  BIFSU  M  PROGRAM 

function  yl=saturl(thl); 

y  1  =Zw*U(n)A2*tan(th  1  )+(W-B  1  )*cos(thl )+Zdlt*U(n)A2*(-0  4) . 

-(0. 5)*rho*CD*  A*U(n)A2*tan(th  1  )*abs(tan(th  1 )); 
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APPENDIX  B 


DTRC  SUBOFF  MODEL  CHARACTERISTICS 
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DESIGN  PARAMETERS  -  SAIL 


Span 

Rooc  Chord 

Tip  Chord 

??  co  Sail  LE 

As pec c  Rac:o 

Discanca 
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0.729  : 
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DESIGN  PARAMETERS   -  CONTROL  PLANE 
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-■ 

CALCULATED   PARAMETER  -  CONTROL  PLANE 


.anzorm  Area 


<£t2> 


0.267 
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Monduneas-Jtial  ozzsezs   and   crass   3ec:.onai   areas    far   :^e 
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